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Time-of-flight  (ToF)  based  acoustic  proximity  ranging  is  widely  used  in  many 
applications.  In  this  dissertation,  its  applications  to  cavity  thickness  monitoring  of 
a supercavitating  vehicle  are  studied.  Most  of  the  currently  available  ToF  based 
acoustic  ranging  systems  are  not  directly  applicable  to  this  case  due  to  their  low 
measurement  accuracy  and  low  parameter  update  rate.  New  measurement  schemes 
and  the  corresponding  signal  processing  approaches  need  to  be  devised  and  their 
performance  evaluated  for  this  challenging  practical  problem.  Based  on  this  motiva- 
tion, four  proximity  ranging  methods,  namely  the  phase-shift  approach,  the  multi- 
frequency technique,  the  PEARS  (Parameter  Estimation  for  Acoustic  Ranging  Sys- 
tems) scheme,  and  the  Multi-PEARS  (Multi-echo  Parameter  Estimation  for  Acoustic 
Ranging  Systems)  algorithm,  are  proposed  and  investigated. 

For  the  phase-shift  approach,  two  frequencies  are  used  and  the  measurements 
are  assumed  to  be  corrupted  by  colored  Gaussian  noise.  By  taking  the  a priori 
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knowledge  of  the  acoustically  hard  reflection  into  account,  a new  time  delay  esti- 
mation algorithm  based  on  the  maximum-likelihood  (ML)  theory  is  derived.  It  is 
shown  that  our  new  method  outperforms  the  traditional  method  in  terms  of  both 
the  estimation  accuracy  and  the  robustness  against  data  model  mismatch.  For  the 
multi-frequency  technique,  an  arbitrary  number  of  frequencies  is  used.  A novel  time 
delay  estimator  based  on  the  nonlinear  least  squares  (NLS)  fitting  criterion  is  derived. 
To  minimize  the  highly  oscillatory  cost  function,  an  efficient  two-stage  estimation  al- 
gorithm is  proposed.  Numerical  examples  show  that  for  a fixed  frequency  interval 
and  a fixed  signal- to-noise  ratio  (SNR),  the  more  frequencies  used,  the  lower  the 
SNR  threshold  and  the  higher  the  estimation  accuracy  that  can  be  obtained.  In- 
spired by  the  multi-frequency  technique,  the  PEARS  scheme  is  devised.  PEARS 
is  novel  in  that  it  is  applicable  to  arbitrary  transmitted  waveforms  as  long  as  the 
waveform  is  periodic.  Numerical  examples  and  the  experiments  performed  by  using 
commercially  available  ultrasonic  transducers  are  used  to  demonstrate  the  excellent 
performance  of  PEARS.  To  deal  with  the  interference  of  secondary  echoes  observed  in 
experiments,  the  Multi-PEARS  algorithm  is  presented  for  the  joint  proximity  rang- 
ing and  secondary  echo  mitigation.  Both  numerical  and  experimental  results  verify 
that  Multi-PEARS  can  provide  very  accurate  distance  measurements  even  in  the 
presence  of  strong  secondary  echoes.  Finally,  ranging  experiments  are  conducted  in 
the  air-water  tunnel  at  different  flow  conditions.  Experimental  results  show  that  the 
distance  between  the  sensors  and  the  air-water  interface  can  be  accurately  measured 
by  using  the  ranging  techniques  developed  in  this  work. 
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CHAPTER  1 
INTRODUCTION 


This  dissertation  is  concerned  with  the  time-of- flight  (ToF)  based  acoustic 
proximity  ranging  approaches  and  their  applications  to  supercavitating  vehicles.  This 
chapter  serves  as  a general  introduction  to  the  background  and  scope  of  the  work. 

1.1  Background 

Measuring  the  distance  between  a known  base  location  and  the  surface  of 
a proximal  object  by  using  an  acoustic  (e.g.,  ultrasound)  signal  is  referred  to  as 
the  acoustic  proximity  ranging.  Most  of  the  existing  acoustic  ranging  systems  use 
the  ToF  based  measurement  method  which  relies  on  the  fact  that  the  propagation 
time  r of  the  ultrasound  is  related  to  the  distance  d from  the  sensor  to  the  object 
by  d — cr/2,  where  c is  the  sound  propagation  velocity  in  the  medium.  Acoustic 
proximity  ranging  is  very  important  in  a wide  range  of  remote  sensing  applications, 
including  nondestructive  testing,  process  monitoring,  robotic  ranging  and  position 
as  well  as  liquid  level  detection  [17,  46,  103].  In  this  dissertation,  we  consider  a 
novel  application  of  acoustic  proximity  ranging  in  monitoring  the  cavity  thickness  of 
a supercavitating  vehicle. 

When  the  flow  of  inviscid  fluid  passes  through  a solid  object  such  as  a cylinder, 
the  flow  speed  around  the  object  is  redistributed.  At  the  regions  where  the  speed 
being  accelerated,  the  corresponding  pressures  on  the  object  surface  decrease.  If  the 
pressure  becomes  less  than  the  vapor  pressure  of  the  fluid,  the  vapor  pockets  are 
formulated.  This  phenomenon  is  referred  to  as  the  natural  cavitation  [9].  Supercavi- 
tation is  an  extreme  version  of  cavitation  characterized  by  the  formation  of  a single 
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bubble  almost  completely  enveloping  the  moving  object  [105].  Recently,  supercavita- 
tion has  found  applications  in  underwater  high  speed  vehicles,  which  are  the  so-called 
supercavitating  vehicles  [4], 

For  an  underwater  vehicle,  the  maximum  speed  is  limited  up  to  about  35  m/s 
by  using  traditional  underwater  techniques.  However,  by  using  the  supercavitation 
technique,  the  speed  of  the  vehicle  can  be  greatly  increased  [4,  51,  57,  93].  Typically, 
at  velocities  over  50  m/s,  the  blunt-nosed  cavitator  and  the  tip-mounted  gas-injection 
system  produce  a low  density  gas  bubble  which  encapsulates  the  entire  body  of  the 
vehicle.  This  bubble  greatly  reduces  the  wetted  friction  drag  and  significantly  in- 
creases the  maximum  speed  of  the  underwater  vehicle  to  around  100  m/s.  Figure 
1.1(a)  shows  a conceptual  drawing  of  such  a supercavitating  vehicle  traveling  at  a 
high  speed.  One  possible  scenario  for  the  supercavitation  is  shown  in  this  figure 
where  the  cavity  interface  surrounding  the  vehicle  can  be  roughly  classified  into  three 
regions.  Near  the  tip  of  the  vehicle,  the  gas-water  interface  has  a smooth,  well-defined 
structure.  In  the  middle  of  the  vehicle,  the  cavity  structure  begins  to  degrade  due 
to  the  instability  of  the  gas-water  interface.  Towards  the  end,  the  process  of  cav- 
ity degradation  continues  downstream  with  the  introduction  of  multiple-bubbles  and 
two-phase  flows. 

The  feasibility  of  such  a vehicle  is  critically  dependent  upon  the  maintenance 
of  the  gaseous  cavity  [53].  As  the  vehicle  travels  in  water  at  high  speeds,  the  cavity 
shape  continuously  changes,  especially  when  the  vehicle  turns.  The  knowledge  of  the 
cavity  thickness  around  the  hull  is  very  important  because  it  provides  the  feedback 
information  to  the  guidance,  control,  and  ventilation  systems  for  maintaining  the 
vehicle  trajectory  in  water.  To  acquire  the  cavity  shape  information,  the  cavity 
sensing  system  needs  to  be  developed.  In  general,  the  cavity  sensing  system  involves 
the  radiation  of  energy  from  the  hull  towards  the  gas- water  interface  and  the  detection 
of  the  reflected  signal.  In  addition,  this  system  should  also  be  1)  sensitive  to  the 
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(a) 


Figure  1.1:  Conceptual  drawing  of  the  cavity  thickness  monitoring  system,  (a)  An 
underwater  supercavitating  vehicle  traveling  at  a high  speed,  (b)  An  acoustic  prox- 
imity ranging  system  mounted  on  the  hull  for  cavity  thickness  monitoring. 

gas-water  interface,  2)  capable  of  measuring  cavity  thickness  of  10  — 100  mm,  3) 
able  to  provide  fast  parameter  updates,  and  4)  small  in  physical  size.  Due  to  the 
acoustically  hard  boundary  at  the  gas- water  interface,  acoustic-based  cavity  thickness 
measurement  methods  are  expected  to  achieve  higher  signal-to-noise  ratio  (SNR)  than 
other  possibilities  such  as  electromagnetic  or  laser-based  methods1.  Meanwhile,  by 
using  micro-electro-mechanical  systems  (MEMS)  techniques,  the  acoustic  transducers 
can  be  fabricated  very  small  in  physical  size  with  the  capability  of  working  in  the 
harsh  seawater  environment  [26].  To  obtain  the  three-dimensional  (3-D)  information 
about  the  gas  cavity,  the  acoustic  cavity  thickness  monitoring  system  is  composed  of 
many  sub-systems  which  are  evenly  mounted  around  the  hull  of  the  supercavitating 
vehicle.  Figure  1.1(b)  shows  an  acoustic  proximity  ranging  sub-system.  The  desired 
goal  of  this  research  is  to  devise  signal  processing  algorithms  for  the  acoustic  proximity 
ranging  system  for  the  real-time  monitoring  of  the  cavity  thickness  of  the  high  speed 
supercavitating  vehicles. 

Generally  speaking,  there  are  two  important  considerations  for  an  acoustic 

proximity  ranging  system.  The  first  consideration  is  the  sensor  properties  including 

1For  the  air- water  interface,  the  power  reflection  coefficients  for  acoustic  signal  and  light  are 
approximate  to  1 and  0.02,  respectively  [14,  47]. 
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resonant  frequency,  bandwidth,  directivity,  efficiency,  and  ruggedness  of  the  ultra- 
sonic transducers  [74,  75].  Specifically,  for  the  working  frequency,  given  a sensing 
range  of  a proximity  ranging  system,  the  working  frequency  should  be  chosen  in  the 
sense  that  the  dissipation  for  that  frequency  is  negligible.  Because  the  dissipation 
of  the  sound  is  proportional  to  the  square  of  the  frequency  [14],  the  larger  the  sens- 
ing range,  the  lower  the  frequency  that  should  be  used.  (For  example,  the  Polaroid 
6500  series  sonar  ranging  system  is  able  to  measure  a distance  up  to  10  m with  the 
working  frequency  of  49.4  kHz  [87].  The  Massa  E-201B/150  has  a ranging  distance 
up  to  1.5  m with  the  working  frequency  of  150  kHz  [73].)  In  addition,  based  on  the 
Rayleigh  criterion  [10],  the  working  frequency  should  also  be  selected  so  that  the 
corresponding  wavelength  is  much  larger  than  the  irregularity  height  on  the  target 
surface.  Otherwise,  only  a very  small  portion  of  the  reflected  signal  can  be  received 
due  to  the  scattering.  The  bandwidth  of  the  transducer  is  an  important  concern  in 
applications  of  range  profile  evaluation  or  multiple  target  detection.  For  these  scenar- 
ios, the  wideband  signals  need  to  be  transmitted  and  received  to  obtain  a high  range 
resolution.  Therefore,  the  transducers  with  a large  bandwidth  are  desired  in  these 
cases  [37,  94,  71].  The  directivity  of  a transducer  crucially  depends  on  ka,  where  k 
is  the  wave  number  and  a is  the  radius  of  the  vibration  piston  or  membrane.  For 
ka  » 1,  the  transducer  can  produce  a beam  pattern  with  a very  strong  directional 
selectivity.  For  b<l,  the  directivity  of  the  transducer  is  omnidirectional  [14]. 

The  second  consideration  is  the  choice  of  signal  waveforms  used  by  the  rang- 
ing system.  Most  of  the  ultrasonic  transducers  can  work  either  in  a pulse  mode  (i.e., 
the  transducer  is  excited  by  a burst  of  pulses),  or  in  a continuous  wave  (CW)  mode 
(i.e.,  the  transducer  is  excited  by  a set  of  continuous  harmonic  signals).  In  the  pulse 
mode,  the  time  delay  (i.e.,  ToF)  is  obtained  by  measuring  the  time  of  arrival  of  the 
reflected  pulses  via  some  very  simple  methods,  such  as  threshold  detection  and  sliding 
window  techniques  [6].  Due  to  its  simplicity,  the  pulse-echo  method  is  widely  used 
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Figure  1.2:  Typical  waveforms  used  by  the  pulse  method  with  the  acoustic  trans- 
ducers. The  waveforms  are  measured  from  Panasonic  ceramic  transducers  with  their 
parameters  listed  in  Table  1.1.  (a)  A sharp  pulse  used  to  excite  the  transmitter,  (b) 
The  response  of  the  receiver  due  to  the  sharp  pulse  in  (a). 

in  commercial  acoustic  ranging  systems.  To  monitor  a fast  changing  environment, 
the  pulse  repetition  interval  should  be  very  small.  However,  most  acoustic  transduc- 
ers are  resonant  devices  [17,  49].  These  kinds  of  acoustic  transducers  are  typically 
underdamped  second-order  dynamic  systems  [96,  36].  Hence,  the  response  of  the 
acoustic  transducer  due  to  a sharp  pulse  inevitably  rings.  The  typical  waveforms 
used  by  the  pulse  method  are  shown  in  Figure  1.2.  The  commercially  available  Pana- 
sonic ceramic  transducers  are  used  in  this  example.  Table  1.1  lists  the  parameters  of 
the  transducers  [86].  Figure  1.2(b)  is  the  transmitted  waveform  measured  when  the 
transmitter  and  the  receiver  are  face  to  face.  It  can  be  noted  from  Figure  1.2  that  the 
pulse  repetition  interval  is  lower  bounded  by  the  length  of  the  pulse  response  (i.e., 
the  transmitted  pulse),  which  limits  its  ability  to  update  the  time  delay  measurement 
quickly.  Moreover,  if  a single  transmitting/receiving  transducer  is  used,  the  problem 
of  the  “dead  zone”  (range  interval  immediately  in  front  of  the  transducer  where  the 
objects  cannot  be  detected)  arises,  which  limits  the  minimum  measurable  distance  of 
the  pulse-echo  method  [71]. 
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Table  1.1:  Parameters  of  the  Panasonic  ultrasonic  ceramic  transducers  used  in  the 
experiments  of  this  work. 


Part  No. 

EFRRSB40K5 

EFRTSB40K5 

Application 

Receiver 

Transmitter 

Resonant  Frequency  (kHz) 

40 

40 

3 dB  Bandwidth  (kHz) 

4 

4 

Sensitivity  (dB)* 

-50 

— 

Sound  Pressure  Level  (dB)** 

— 

105 

Maximum  Input  Voltage  (Vrms) 

— 

20 

* 0 dB=l  V/Pa  **  0 dB=2  x 1(T5  Pa 


In  the  CW  mode,  the  ranging  system  is  designed  to  be  simultaneously  trans- 
mitting and  receiving  signals  by  a transmitter  and  a receiver  separately  [61,  58,  112, 
77],  which  can  be  realized  by  1)  using  several  transducer  pairs  with  different  resonant 
frequencies,  2)  sweeping  the  resonant  frequency  of  one  transmitter  [49],  or  3)  excit- 
ing one  transmitter  with  a set  of  sinusoidal  signals  simultaneously  or  sequentially 
in  a certain  frequency  range.  The  time  delay  information  is  obtained  by  measuring 
the  phase-shift  at  each  frequency  in  the  received  signal  thereafter.  This  measure- 
ment technique  theoretically  does  not  have  “dead  zone”  problems  2 and  can  be  used 
to  obtain  highly  accurate  measurements  even  with  a narrowband  transducer  [61]. 
Moreover,  by  using  the  CW  mode,  the  peak-to-average  ratio  of  the  transmitted  sig- 
nal can  be  effectively  reduced  as  well.  Although  special  cases  of  using  two  or  three 
different  frequencies  have  been  considered  before  [58,  61,  112],  most  of  the  previous 
approaches  do  not  consider  a priori  information  of  the  reflection  coefficient  for  the 
reflection  boundary.  Hence,  the  corresponding  time  delay  estimation  methods  are 
not  optimal  in  terms  of  the  estimation  accuracy.  In  addition,  no  efficient  ways  were 
presented  before  to  measure  the  phase-shift  at  each  frequency  with  the  consideration 
for  the  real-time  implementation. 

2In  practice,  the  minimum  detectable  distance  for  the  CW  system  can  be  limited  by  other  issues 
such  as  the  directivities  of  the  transducers. 
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In  conclusion,  the  new  proximity  range  measurement  schemes  and  the  corre- 
sponding signal  processing  approaches  need  to  be  developed  for  the  real-time  moni- 
toring of  the  cavity  thickness.  This  is  the  motivation  for  our  work. 

1.2  Scope  of  the  Work 

Compared  with  the  pulse  mode,  the  ranging  system  working  in  the  CW  mode 
does  not  have  “dead  zone”  problems  and  can  be  used  to  obtain  high  measurement 
accuracy  by  using  the  narrowband  transducer.  Therefore,  the  CW-based  distance 
measurement  method  is  expected  to  be  superior  to  its  pulse-based  counterpart.  Due 
to  the  above  reasons,  this  work  starts  with  the  CW-based  method.  The  simplest  CW- 
based  method  is  the  so-called  phase-shift  approach  in  which  only  two  frequencies  are 
used.  The  phase-shift  approach  is  widely  used  in  optics  [15,  52,  91].  A traditional 
estimator  for  the  phase-shift  approach  is  described  in  Hornung  and  Brand  [49].  It  is 
based  on  the  fact  that  given  two  different  working  frequencies,  the  difference  between 
the  two  phase-shifts  is  proportional  to  the  time  delay.  However,  this  estimator  is  not 
optimal  when  the  signal  is  reflected  from  the  boundary  with  a known  reflection  coef- 
ficient. Furthermore,  this  estimator  does  not  take  the  noise  into  account.  In  Chapter 
3 , the  phase-shift  approach  is  studied.  Based  on  the  assumptions  that  the  amplitude 
of  the  signal  is  a nonnegative  unknown  and  the  measurements  are  corrupted  by  an 
unknown  colored  Gaussian  noise,  a new  time  delay  estimator  based  on  the  maximum- 
likelihood  (ML)  approach  is  derived.  However,  this  estimator  requires  the  solution 
of  an  optimization  problem  whose  objective  function  is  highly  oscillatory.  Inspired 
by  related  work  [110],  we  propose  an  approach  to  obtain  an  initial  time  delay  esti- 
mate by  assuming  that  the  amplitude  of  the  reflected  signal  is  complex-valued  and 
then  to  refine  the  initial  estimate  based  on  the  true  cost  function.  As  shown  in  this 
work,  the  root-mean-squared  errors  (RMSEs)  of  the  parameter  estimates  of  the  new 
method  can  approach  the  corresponding  Cramer- Rao  bounds  (CRBs),  which  are  the 
best  performance  bounds  for  any  unbiased  estimators,  as  the  SNR  increases  [24,  99]. 
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In  Chapter  4,  the  discussion  is  extended  to  a multi-frequency  case,  in  which 
an  arbitrary  number  of  frequencies  can  be  used.  However,  the  most  difficult  problem 
in  the  practical  implementation  is  to  measure  the  phase-shifts  for  many  frequencies 
efficiently.  As  will  be  shown  in  Chapter  4,  this  difficulty  is  overcome  by  our  specially 
designed  excitation  waveform  in  which  the  signals  with  different  frequencies  are  or- 
thogonal to  each  other  for  a given  time  period.  Then,  we  formulate  a general  data 
model  applicable  to  this  multi-frequency  measurement  scheme  where  the  amplitude 
of  the  reflected  signal  is  considered  to  be  a nonnegative  unknown.  According  to  this 
data  model,  a new  time  delay  estimator  based  on  the  nonlinear  least  squares  (NLS) 
fitting  criterion  is  derived.  To  optimize  the  NLS  cost  function,  a two-stage  fast  op- 
timization algorithm  is  presented.  Numerical  results  show  that  this  method  avoids 
the  search  over  the  parameter  space,  yet  there  are  no  obvious  performance  degrada- 
tions when  compared  to  the  direct  search  method.  Numerical  results  also  show  that 
the  RMSEs  of  the  new  method  can  approach  the  corresponding  CRBs  as  the  SNR 
increases. 

In  Chapter  5,  inspired  by  the  parameter  estimation  approach  adopted  by  the 
multi-frequency  technique,  the  PEARS  (Parameter  Estimation  for  Acoustic  Ranging 
Systems)  scheme  is  devised.  It  will  be  shown  that,  after  applying  the  fast  Fourier 
transform  (FFT)  to  one  period  of  the  sampled  signal  with  an  arbitrary  waveform,  in 
the  frequency  domain,  a data  model  similar  to  that  used  by  the  multi-frequency  tech- 
nique can  be  formulated.  Hence,  the  two-stage  fast  time  delay  estimation  algorithm 
developed  for  the  multi-frequency  technique  can  be  applied  to  this  new  case  with  only 
minor  modifications.  The  ability  to  deal  with  the  time  delay  estimation  problem  for 
an  arbitrary  transmitted  waveform  provides  more  flexibility  in  designing  a practical 
proximity  ranging  system.  Then,  a numerical  example  is  used  to  illustrate  the  perfor- 
mance of  PEARS.  Numerical  results  show  that  the  RMSEs  of  the  estimates  obtained 
by  PEARS  can  achieve  the  corresponding  CRBs  as  the  SNR  increases.  The  efficiency 
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of  PEARS  is  also  demonstrated  by  performing  experiments  using  the  commercially 
available  transducers. 

One  problem  noticed  in  experiments  is  the  interference  of  the  secondary  echoes. 
When  ultrasonic  transducers  face  a sound-hard  reflective  surface,  the  reflected  sound 
wave  can  bounce  back  and  forth  several  times  between  the  sensors  and  the  reflection 
surface  before  decaying  to  zero,  which  results  in  unwanted  strong  and  overlapping  sec- 
ondary echoes  in  the  received  signal.  The  presence  of  the  strong  secondary  echoes  can 
degrade  the  performance  of  our  former  developed  methods.  To  solve  this  problem,  in 
Chapter  6,  the  Multi-PEARS  (Multi-echo  Parameter  Estimation  for  Acoustic  Rang- 
ing Systems)  algorithm  is  presented  for  the  joint  proximity  ranging  and  secondary 
echo  mitigation.  Both  numerical  and  experimental  results  show  that  Multi-PEARS 
can  provide  very  accurate  distance  measurements  even  in  the  presence  of  strong  sec- 
ondary echoes. 

Finally,  in  Chapter  7,  more  realistic  experiments  are  performed  in  the  air- 
water  tunnel  under  different  flow  conditions.  Three  different  air-water  interfaces 
(i.e.,  stationary,  nonstationary,  and  nonstationary  with  periodic  disturbances)  are 
measured.  Experimental  results  show  that  the  distance  between  the  sensors  and 
the  air-water  interface  can  be  accurately  measured  by  the  novel  acoustic  proximity 
ranging  techniques  developed  in  this  work. 

1.3  Organization  of  the  Dissertation 

This  dissertation  is  organized  as  follows.  Chapter  2 gives  a literature  survey 
of  such  topics  as  proximity  sensing,  acoustic  proximity  ranging,  and  time  delay  esti- 
mation. In  Chapter  3,  a new  time  delay  estimator  for  the  phase-shift  approach  based 
on  the  ML  theory  is  derived.  The  model  mismatch  analysis  of  this  new  time  delay 
estimator  and  the  numerical  examples  are  also  presented.  Chapter  4 investigates  the 
multi-frequency  technique,  in  which  an  arbitrary  number  of  excitation  frequencies 
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can  be  used.  In  Chapter  5,  the  PEARS  scheme  is  devised.  Both  numerical  and  ex- 
perimental results  are  used  to  demonstrate  the  performance  of  PEARS.  In  Chapter  6, 
the  Multi-PEARS  algorithm  is  presented  to  deal  with  the  secondary  echoes.  Chapter 

7 shows  the  experimental  results  obtained  in  the  air-water  tunnel.  Finally,  Chapter 

8 concludes  this  dissertation. 


CHAPTER  2 
LITURATURE  SURVEY 


One  of  the  most  important  techniques  for  proximity  sensing  is  the  acoustic 
proximity  ranging.  This  field  of  research  has  been  well  documented  in  the  literature. 
Many  excellent  articles,  both  theoretical  and  application-oriented,  are  available  (for 
example,  [48,  61,  112,  77,  108,  68,  56,  6,  44,  16,  79]  and  the  references  therein). 
Classical  works  and  typical  application  examples  have  been  edited  and  reprinted  as 
well  [46].  This  chapter  gives  a brief  review  of  the  subjects  that  are  related  to  our  work, 
including  proximity  sensing,  acoustic  proximity  ranging,  and  time  delay  estimation. 

2.1  Proximity  Sensing 

The  purpose  of  this  section  is  to  provide  an  overview  of  proximity  sensing 
techniques.  This  section  helps  us  understand  where  the  acoustic  proximity  ranging 
lies  among  all  the  proximity  sensing  techniques.  In  addition,  it  can  also  help  us  un- 
derstand why  or  why  not,  in  some  practical  applications,  acoustic  proximity  ranging 
is  a better  choice  than  the  others. 

Proximity  sensing  can  be  broken  down  into  two  categories  including  the  pres- 
ence detection  and  the  distance  measurement.  The  most  often  used  techniques  are 
illustrated  in  Figure  2.1.  Many  of  them  are  commercially  available  as  well  [41]. 

2.1.1  Presence  Detection 

For  the  presence  detection,  the  sensing  systems  are  primarily  designed  to 
determine  the  presence  of  a proximal  object.  Field-based  [36,  20,  55,  31,  92,  89,  30], 
photoelectric  [36,  38,  20],  and  acoustic  sensing  [18,  59,  50,  20]  systems  are  commonly 
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Figure  2.1:  Proximity  sensing  techniques. 


used  for  this  purpose.  We  concentrate  on  the  distance  measurement  approaches  due 
to  the  aforementioned  supercavitation  applications. 

2.1.2  Distance  Measurement 

For  the  distance  measurement,  the  sensing  systems  are  designed  not  only 
to  detect  the  object,  but  also  to  provide  the  distance  between  the  sensor  and  the 
object.  The  distance  measurement  methods  can  be  cast  into  three  categories  as 
signal  intensity,  triangulation,  and  ToF  based  methods. 

Signal  intensity.  Two  signal  intensity  based  distance  ranging  methods 
have  been  reported  before.  The  first  method  is  the  so-called  amplitude  method  [49], 
which  can  be  analyzed  in  a way  similar  to  the  Fabry-Perot  principle  known  in  optics 
[47].  In  the  space  between  the  sensor  and  object,  the  emitted  and  reflected  waves 
interfere  with  each  other.  The  received  signal  amplitude  is  maximal  or  minimal  when 
waves  are  in  resonance  or  antiresonance.  The  distance  between  the  sensor  and  the 
object  is  obtained  by  measuring  the  frequency  of  occurrence  between  two  neighboring 
resonance  peaks  or  notches.  The  second  method  is  based  on  the  inverse  square  law  of 
the  emitted  energy,  which  states  that  as  the  distance  from  a point  source  increases, 
the  intensity  of  the  source  diminishes  as  a function  of  the  square  of  the  distance.  To 
remove  the  uncertainty  of  the  intensity  of  the  reflected  signal  due  to  the  different 
target  surfaces,  a scheme  described  in  Nansel  [84]  makes  use  of  a pair  of  identical 
point-source  light-emitting  diodes  (LEDs)  positioned  at  a known  distance  apart,  with 
their  incident  light  focused  on  the  target  surface.  A photodetector  is  placed  in  the 
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middle  of  the  two  LEDs  as  a receiver.  The  difference  between  the  received  signal 
intensities  corresponding  to  the  different  light  sources  can  be  exploited  to  yield  the 
absolute  range  value.  The  advantage  of  the  signal  intensity  based  methods  lies  in 
their  simplicity  in  both  the  concept  and  implementation.  However,  their  use  is  limited 
due  to  the  low  measurement  accuracy  compared  with  the  triangulation  and  the  ToF 
based  methods  discussed  below.  In  addition,  it  should  be  noted  that,  in  acoustics,  the 
inverse  square  law  is  only  valid  when  the  distance  from  the  transducer  is  greater  than 
ka  [14],  where  k is  the  wave  number  and  a is  the  radius  of  the  piston  or  membrane 
in  the  transducer. 

Triangulation.  Triangulation  presents  a simple  trigonometric  method  for 
calculating  the  distance  and  angles  needed  to  determine  the  object  location  [36,  20]. 
The  triangulation  ranging  method  is  based  on  the  fact  that  given  the  length  of  one 
side  and  two  angles  of  a triangle,  it  is  possible  to  determine  the  lengths  of  the  other 
two  sides  and  the  third  angle.  Nowadays,  the  triangulation  method  is  widely  used  in 
optical  displacement  sensors  and  can  provide  very  accurate  distance  measurements 
[21,  72,  29].  However,  a major  shortcoming  of  the  triangulation  principle  based 
optical  displacement  sensors  is  that  they  usually  have  a relatively  small  measurement 
range  and  require  a long  stand  off  distance  [21],  For  example,  the  ILD  2000-10 
laser  displacement  sensor  [78]  (which  will  be  used  in  our  experiment  later  on)  has  an 
accuracy  of  0.5  gm  with  the  measurement  range  of  10  mm  and  the  stand  off  distance 
of  58  mm. 

Time-of-flight.  The  ToF  based  ranging  systems  measure  the  round  trip 
travel  time  between  the  transmitted  and  reflected  signal  from  the  object.  The  distance 
between  the  sensor  and  the  object  is  determined  by  d,  = cr/2,  where  c is  the  velocity 
of  the  energy  wave  in  the  medium  and  r is  ToF.  In  this  section,  we  address  only  the 
laser  and  the  microwave  based  systems.  The  ultrasonic  ToF  ranging  system  will  be 
addressed  in  detail  in  the  next  section. 
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The  ToF  ranging  system  using  laser  is  referred  to  as  lidar  (light  detection 
and  ranging).  Due  to  the  narrow  beamwidth  of  the  laser  light,  lidar  is  capable  of 
producing  a very  accurate  object  resolution  [107,  1,  98,  60].  The  only  problem  with 
lidar  is  the  difficulty  in  measuring  the  very  small  ToF.  For  general  pulse-echo  lidar 
systems,  an  accuracy  of  1.5  m requires  the  time  measurement  of  one  nanosecond, 
which  needs  rather  sophisticated  electronic  detection  circuitry  [107].  One  solution 
to  this  problem  is  to  measure  the  phase-shifts  of  a continuous  light  wave  or  a long 
modulated  light  pulse  rather  than  directly  measure  the  ToF  of  a short  light  pulse 
[15,  52,  91]. 

A good  example  of  the  proximity  ranging  system  by  using  microwave  is  the 
micro  power  impulse  radar  (MIR)  developed  by  the  Lawrence  Livermore  National 
Laboratory  [5].  The  fundamental  operating  principle  of  MIR  is  the  same  as  a con- 
ventional pulsed  radar.  However,  MIR  makes  use  of  a noise  generator  to  trigger  a 
pulse  generator,  which  can  produce  very  short  and  very  low-power  pulses.  A very 
high-speed  sampling  receiver  is  used  to  sample  the  received  pulses.  The  sampled 
signals  are  passed  through  an  averaging  electronic  circuit,  which  averages  10,000 
received  pulses  to  obtain  one  ToF  measurement.  The  measurement  range  for  the 
standard  MIR  is  about  6 m.  The  most  striking  features  of  the  MIR  system  are  its 
low  cost  and  low  power  consumption  [36].  Compared  with  the  triangulation  based 
ranging  systems,  the  ToF  based  ranging  systems  provide  slightly  lower  measurement 
accuracy.  However,  their  measurement  range  is  much  larger. 

A brief  summary  of  the  different  distance  measurement  methods  is  listed  in 
Table  2.1. 

2.2  Acoustic  Proximity  Ranging 

In  this  section,  a brief  overview  of  the  properties  of  ultrasonic  transducers  is 
presented,  followed  by  a review  of  the  various  ToF  measurement  schemes  used  by 
commercial  and  experimental  acoustic  ranging  systems. 
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Table  2.1:  Summary  of  the  different  distance  measurement  methods. 


Methods 

Principles 

Merits 

Limitations 

Signal 

Intensity 

Amplitude 

Method 

Fabry-Perot  principle. 

Simplicity;  Widely 
used  in  optics. 

Need  to  sweep  a 
certain  frequency 
range. 

Inverse 
Square  Law 

Intensity  of  received 
signal  deceases  as 
function  of  ct . 

Easy  to  implement. 

Sensitive  to 
reflection  surface; 
Low  accuracy;  In 
acoustics,  d » ka . 

Triangulation 

Trigonometric 

relationship. 

Very  accurate; 
Widely  used  in 
optics. 

Expensive; 
Relatively  small 
measurement  range; 
Requires  standoff 
distance. 

ToF 

Lidar,  Radar 

d = ct  12,  where  cis 
EM/sound  wave 
speed  and  r is  time 
delay. 

High  accuracy; 
Long  ranging 
distance. 

Rather 

sophisticated 

systems. 

Sonar 

Rather  simpler 
systems;  Robustness 
against  EM  noise  and 
interference. 

Need  to  know  c. 

2.2.1  Ultrasonic  Transducers 

Transducers  are  very  important  components  for  any  acoustic  ranging  systems. 
Their  temporal  and  spatial  characteristics  are  prime  determinants  to  the  performance 
of  a ranging  system.  A detailed  survey  of  ultrasonic  transducers  for  use  in  air  is  pro- 
vided by  Manthey  et  al.  [71].  Various  kinds  of  commercial  ultrasonic  transducers  for 
different  application  purposes  are  also  available  (for  example,  [73,  87,  79]).  Ultra- 
sonic transducers  are  categorized  by  the  physical  principles  utilized  to  generate  the 
ultrasound,  such  as  electrostatic,  electrodynamic,  piezoelectric,  and  magnetrostric- 
tive  principles.  For  industrial  applications,  electrostatic  and  piezoelectric  transducers 
are  of  main  interest. 

Electrostatic  transducers  are  similar  to  electrical  capacitors  which  consist  of 
one  fixed  plate  and  one  movable  plate.  A thin  metal  sheet  or  a single-sided  metallized 
plastic  foil  is  used  as  a diaphragm  to  convert  electrical  energy  into  acoustic  energy 
and  vice  versa  [102,  39,  43,  85,  37,  94].  The  advantages  of  electrostatic  transducers 
include  high  transfer  coefficients  and  wide  bandwidths  in  both  the  transmitting  and 
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(a)  (b) 


Figure  2.2:  Conventional  and  UF  MEMS-based  ultrasonic  transducers,  (a)  Photo- 
graph of  three  examples  of  conventional  transducers,  (b)  Top  view  of  the  microscopic 
image  of  the  MEMS-based  transducer  developed  by  IMG  at  UF  [26]. 

the  receiving  modes.  However,  their  use  is  limited  because  of  the  sensitivity  to  dust 
and  humidity,  low  mechanical  ruggedness,  and  the  need  for  an  external  bias  voltage 
[68]. 

Clearly,  piezoceramic  transducers  dominate  the  modern  industrial  ultrasonic 
measurement  systems  [36,  43].  As  the  piezoelectric  material  is  excited  by  an  electrical 
potential,  the  physical  size  of  the  piezoelectric  material  expands  or  contracts.  The 
moving  cone  displaces  air  and  generates  ultrasound.  Because  the  piezoelectricity  is 
a reversible  phenomenon,  the  piezoelectric  material  generates  the  voltage  as  the  in- 
coming ultrasonic  wave  flexes  it.  Thus,  ultrasound  can  be  detected  by  using  the  same 
device  as  well.  Although  piezoelectric  transducers  are  usually  less  sensitive  than  elec- 
trostatic transducers,  the  typical  advantages  are  their  compact,  rugged  mechanical 
design,  good  impedance  matching,  and  the  wide  range  of  the  operation  temperature 
[68]. 

Conventional  transducers  based  on  the  aforementioned  ultrasound  generating 
approaches  typically  have  a large  physical  size,  with  a transducer  diameter  ranging 
from  1 cm  to  5 cm.  If  the  driving  and  read-out  circuits  are  taken  into  account,  the 
housing  geometry  for  the  entire  system  usually  exceeds  6x3x2  cm3  [49].  Figure  2.2(a) 
shows  a photograph  of  three  types  of  conventional  ultrasound  proximity  transducers. 
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MEMS-based  transducers  have  been  developed  using  a combination  of  stan- 
dard semi-conductor  processing  and  silicon  micromachining  techniques  [96].  These 
transducers  can  offer  several  potential  advantages  over  their  conventional  counter- 
parts. The  MEMS-based  transducers  can  be  fabricated  in  a very  small  size  (on  the 
order  of  1 mm)  [32,  83].  By  using  the  batch  fabrication  process,  transducers  can 
be  made  with  better  matched  characteristics  and  much  lower  cost  per  device,  which 
would  be  very  suitable  for  a distributed  sensor  network  [102,  3,  85].  Moreover,  the 
integrated  components  reduce  the  size  of  the  housing  and  increase  the  robustness 
of  the  entire  system.  Micromachined  ultrasonic  transducers  based  on  different  op- 
eration principles  have  been  reported  including  electrostatic,  thermal-piezoresistive 
and  piezoelectric  principles  [81,  45].  In  our  particular  application,  electro-thermal 
actuation  and  piezoresistive  sensing  are  adopted  because  they  provide  the  best  com- 
promise between  performance,  cost,  durability  and  ease  of  fabrication  [49].  The  first 
generation  of  the  MEMS-based  acoustic  proximity  transducer  developed  by  the  inter- 
disciplinary microsystems  group  (IMG)  at  the  University  of  Florida  (UF)  is  shown  in 
Figure  2.2(b)  [26].  The  diaphragm  of  this  transducer  is  1 mm  in  diameter  and  10  /rm 
in  thickness,  which  incorporates  a central  resistive  heater  for  thermoelastic  actuation 
and  diffused  piezoresistors  for  sensing  acoustic  pressure. 

Most  of  the  conventional  and  MEMS-based  ultrasonic  transducers  used  in  air 
are  resonant  devices.  These  kinds  of  transducers  are  typically  underdamped  second- 
order  dynamic  systems  (i.e.,  bandpass  systems),  and  hence,  the  response  of  the  acous- 
tic transducer  due  to  the  shape  pulse  inevitably  rings  [8,  6].  The  bandpass  nature 
of  the  ultrasonic  transducers  and  the  corresponding  transmitted  signals  will  be  the 
main  concern  for  our  signal  processing  algorithms. 

2.2.2  ToF  Measurement  Schemes 

In  this  section,  we  briefly  review  the  ToF  measurement  techniques  adopted 
by  commercial  and  experimental  acoustic  ranging  systems.  Most  of  these  methods 
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are  studied  mainly  from  the  implementation  and  application  oriented  standpoints. 
Hence,  they  may  not  be  optimal  in  the  theoretical  sense.  The  ToF  measurement 
problem  will  be  further  reviewed  from  the  statistical  parameter  estimation  theory 
point  of  view  in  the  next  section. 

Pulse-echo  method.  Due  to  its  simplicity  and  reliability,  the  pulse-echo 
together  with  the  threshold  detection  method  is  the  most  widely  used  ToF  mea- 
surement method  in  the  current  commercial  acoustic  ranging  systems.  Two  typical 
examples  are  the  Polaroid  6500  series  sonar  ranging  system  [87]  and  the  Massa  E- 
201B/150  [73]  ranging  system.  For  Polaroid  6500,  a single  electrostatic  transducer  is 
used  for  both  the  transmitter  and  the  receiver.  The  excitation  pulse  consists  of  16 
cycles  at  a frequency  of  49.4  kHz.  To  avoid  the  ringing  tail  of  the  current  transmit- 
ted pulse  being  detected  as  a returned  signal,  an  internal  blanketed  circuit  is  used  to 
block  the  detection  for  a short  period  right  after  the  transmission.  It  is  claimed  that 
Polaroid  6500  is  able  to  measure  a distance  from  0.15  m to  10  m with  an  accuracy 
of  ±1%  of  the  entire  range.  A notable  feature  of  this  system  is  its  multiple  echo 
working  mode,  in  which  the  system  can  resolve  the  echoes  for  objects  placed  0.076 
m apart.  For  Massa  E-201/150,  two  transducers  working  at  150  kHz  are  used  to 
serve  as  the  transmitter  and  the  receiver  separately.  This  bistatic  configuration  can 
provide  a shorter  “dead  zone”  of  0.12  m rather  than  0.20  m for  its  single  transducer 
counterpart.  It  is  claimed  that  the  operating  range  for  this  system  is  from  0.12  m to 
1.5  m with  an  accuracy  of  ±1  mm. 

However,  as  pointed  out  by  Barshan  [6],  the  main  drawback  of  the  threshold 
detection  method  is  that,  on  average,  the  ToF  measurement  is  larger  than  the  actual 
ToF,  which  is  a consequence  of  the  relatively  long  rise  time  of  the  pulse  response 
for  the  low  bandwidth  ultrasonic  transducer.  Barshan  and  Ayrulu  [6,  7],  studied 
four  different  range  measurement  methods  for  pulse-echo  ranging  systems  includ- 
ing threshold,  curve  fitting,  sliding  window,  and  correlation.  For  the  curve  fitting 
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method,  a parabolic  curve  is  used  to  fit  the  signal  envelope  around  the  rising  edge 
of  the  echo.  Based  on  the  vertex  of  the  parabola,  a time  delay  estimate  can  be  ob- 
tained. For  the  sliding  window  method,  a window  of  width  N is  slid  through  the 
echo  signal  one  sample  at  a time.  At  each  window  position,  the  number  of  samples 
which  exceed  a threshold  is  counted.  If  this  number  exceeds  a predetermined  value, 
then  a target  is  assumed  to  be  presented  and  a ToF  estimate  is  produced.  Simulation 
and  experimental  results  show  that  the  correlation  method  is  the  best  in  terms  of 
the  accuracy.  The  curve  fitting,  sliding  window,  and  threshold  methods  follow  the 
correlation  method  in  the  order  of  decreasing  complexity  and  can  offer  an  acceptable 
performance  at  a much  lower  computational  cost. 

CW-based  method.  The  CW-based  method  is  another  approach  com- 
monly used  to  measure  ToF.  In  the  CW-based  method,  two  transducers,  which  serve 
as  the  transmitter  and  the  receiver  separately,  are  used  to  continuously  transmit  and 
receive  the  ultrasound  signal.  The  advantage  of  the  CW-based  method  is  that  there 
is  no  “dead  zone”  problem  as  in  the  pulse-echo  method.  However,  the  crosstalk  may 
be  induced  when  the  transmitter  and  the  receiver  are  closely  placed1,  and  thus  the 
isolation  between  the  transducers  becomes  a key  issue  in  the  real  implementation. 

One  simple  CW-based  ranging  scheme  is  to  transmit  two  or  three  sinusoidal 
signals  at  different  frequencies.  The  phase-shifts  between  the  received  signal  and  the 
transmitted  signal  at  each  frequency  are  measured,  and  ToF  is  obtained  based  on 
the  phase-shift  measurements.  Lee  et  al.  [61]  studied  a CW-based  acoustic  ranging 
system  by  using  three  different  frequencies,  which  are  modulated  around  the  resonant 
frequency  of  the  transducer.  It  is  claimed  that  an  accuracy  of  ±0.3  mm  is  obtained 
over  a 500  mm  range  by  using  the  conventional  narrowband  transducers.  Based  on 
this  framework,  the  phase  errors  and  their  effects  on  the  ToF  estimation  are  discussed 

in  Yang  et  al.  [112].  A two-frequency  CW  range  finder  implemented  by  using  MEMS 

1In  real  systems,  the  crosstalk  can  be  induced  by  many  facts  such  as  the  system  package  and  the 
directivities  of  the  transducers. 
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transducers  is  reported  in  Kuratli  and  Huang  [58].  An  accuracy  of  4 mm  is  obtained 
for  a distance  from  10  mm  to  100  mm.  A similar  experimental  work  by  using  three 
frequencies  with  MEMS  transducers  is  reported  in  Hornung  and  Brand  [49].  However, 
none  of  the  above  methods  present  efficient  ways  to  measure  the  phase-shifts  with 
the  consideration  for  the  real-time  implementation.  In  addition,  they  do  not  take  the 
a priori  information  of  the  acoustically  hard  reflection  into  account. 

The  linear  frequency  modulation  scheme  [76,  77]  is  another  CW-based  ranging 
approach,  which  is  based  on  a continuous  transmission  of  a signal  whose  frequency  is 
linearly  modulated  over  a certain  frequency  band.  The  received  signal  is  mixed  with 
the  transmitted  signal,  and  the  frequency  of  the  output  signal  is  proportional  to  the 
distance  of  the  object.  By  applying  the  Fourier  transform  to  the  output  signal,  the 
range  profile  of  the  object  can  be  obtained.  However,  this  method  requires  special 
transducers  that  can  work  over  a wide  frequency  band.  In  addition,  the  nonlinearity 
of  the  frequency  modulation  induced  by  many  factors  such  as  the  transfer  functions 
of  the  transducers  needs  to  be  compensated  out  accurately  [90]. 

2.3  Time  Delay  Estimation 

The  ToF  measurement,  in  essence,  is  the  time  delay  estimation  problem.  In- 
deed, time  delay  estimation  is  a well-known  problem  that  arises  frequently  in  radar, 
sonar,  radio  navigation,  and  medical  imaging  [22]. 

In  the  single  path  and  white  Gaussian  noise  case,  the  matched  filter  approach 
is  the  optimal  method  [22,  23].  The  matched  filter  method  can  be  implemented 
in  either  the  time  or  the  frequency  domain.  The  frequency  domain  is  preferred  if 
FFT  can  be  used  and  the  transmitted  signal  is  band-limited.  This  is  because  the 
time  domain  correlator  needs  interpolation,  e.g.,  parabolic  fitting  [80],  to  obtain  the 
subsample  resolution.  However,  this  inconvenience  can  be  avoided  by  transforming 
the  data  model  into  the  frequency  domain  where  the  unknown  time  delay  can  take  on 
a continuum  of  values  [66].  In  the  presence  of  multipath,  the  ability  of  the  matched 
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filter  approach  to  estimate  the  time  delays  can  be  affected  by  closely  spaced  echoes. 
This  is  due  to  the  fact  that  the  matched  filter  approach  cannot  resolve  two  signals 
with  a time  spacing  that  is  less  than  the  reciprocal  of  the  signal  bandwidth2  [33]. 

Recently,  many  super-resolution  time  delay  estimation  techniques  have  been 
devised.  In  the  time  domain,  the  received  signal  can  be  modeled  as  the  sum  of 
multiple  scaled  and  delayed  replicas  of  the  transmitted  signal  plus  noise.  In  the 
frequency  domain,  this  data  model  becomes  the  sum  of  multiple  weighted  complex 
exponentials  plus  noise.  The  frequency  domain  model  is  similar  to  those  used  in 
the  sinusoidal  parameter  and  angle  estimation  problems  except  that  the  complex 
exponentials  are  weighted  by  the  known  signal  spectrum.  Based  on  this  observation, 
many  existing  sinusoidal  frequency  and  angle  estimation  algorithms,  such  as  the 
multiple  signal  classification  (MUSIC)  [95],  linear  prediction  [104],  and  ML  methods, 
are  suggested  to  solve  this  problem  [19,  12,  54].  However,  they  are  best  suited  for 
the  complex- valued  signals  with  special  shapes  (such  as  flat  band-limited  spectrum). 
A computationally  efficient  approach  based  on  the  expectation  maximization  (EM) 
algorithm  [82]  is  proposed  in  Feder  and  Weinstein  [35]  that  decouples  a complicated 
multidimensional  optimization  problem  into  a sequence  of  multiple  separate  one- 
dimensional (1-D)  optimization  problems.  However,  the  EM  method  is  very  sensitive 
to  initial  conditions  and  no  systematic  initialization  method  is  given.  A special  case 
of  the  time  delay  estimation  is  when  the  signal  is  bandpass  with  a real-valued  gain, 
which  is  usually  encountered  in  active  sonar  systems  [27].  In  this  case,  the  correlation 
function  between  the  received  signal  and  the  transmitted  signal  will  oscillate  near  the 
carrier  frequency  of  the  transmitted  signal.  As  a consequence,  many  existing  time 
delay  estimation  algorithms  perform  poorly  due  to  the  problem  of  converging  to  local 

optimum  points.  This  challenging  issue  is  addressed  in  Bell  et  al.  [11,  70,  13],  where 

2 For  example,  by  using  the  Panasonic  ultrasonic  transducers  given  in  Table  1.1,  the  minimum 
resolvable  time  spacing  between  two  echoes  is  (4  kHz)-1  = 0.25  ms. 
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the  algorithms  proposed  are  all  based  on  the  NLS  fitting  criterion  and  differ  from 
each  other  in  the  way  the  NLS  cost  function  is  optimized. 

A weighted  Fourier  transform  and  RELAXation  based  method  (WRELAX)  is 
proposed  by  Li  and  Wu  in  [66]  to  minimize  a frequency  domain  NLS  cost  function. 
The  most  striking  feature  of  the  WRELAX  algorithm  is  that  it  decouples  the  mul- 
tidimensional optimization  problem  into  a series  of  1-D  optimization  problems  in  a 
conceptually  and  computationally  simple  way.  The  resolution  of  WRELAX  is  much 
higher  than  the  conventional  matched  filter  approach.  WRELAX  is  also  applied  to 
deal  with  the  problem  that  the  cost  function  is  highly  oscillatory  [110],  in  which  the 
oscillation  problem  is  avoided  by  modeling  the  amplitudes  of  each  received  signal  as 
complex-valued  to  obtain  a good  initial  estimate.  This  is  similar  to  the  methodology 
adopted  in  Vaccaro  et  al.  [106,  70].  To  increase  the  convergence  speed  of  WRELAX 
when  signals  are  very  closely  spaced  in  the  arrival  times,  a method  which  combines 
MODE  (Method  of  Direction  Estimation  [101,  100])  with  WRELAX  is  proposed  in 
Wu  et  al.  [111].  MODE-WRELAX  can  be  used  for  both  complex-  and  real- valued 
signals.  It  outperforms  MODE  in  terms  of  robustness  against  data  model  errors  and 
provides  a better  resolution  than  WRELAX. 


CHAPTER  3 

PHASE-SHIFT  APPROACH 

3.1  Introduction 

CW-based  methods  require  that  the  transducer  generate  ultrasound  continu- 
ously and  receive  the  reflected  signal  simultaneously.  One  of  the  widely  used  CW- 
based  methods  is  referred  to  as  the  phase-shift  approach , which  is  simple  and  fast 
since  it  requires  only  two  measurements  of  a phase-shift  at  two  different  working  fre- 
quencies. A traditional  estimator  for  the  phase-shift  approach  is  described  in  [49]  and 
is  based  on  the  fact  that  given  the  two  different  working  frequencies,  the  difference 
between  the  two  phase-shifts  is  proportional  to  the  time  delay.  However,  it  turns  out 
that  this  estimator  is  not  optimal  when  the  signal  is  known  to  have  a real-valued 
or  nonnegative  amplitude.  Furthermore,  this  estimator  does  not  take  the  noise  into 
account. 

In  this  chapter,  we  consider  the  phase-shift  approach  in  the  case  when  the  am- 
plitude of  the  signal  is  a nonnegative  unknown,  and  the  measurements  are  corrupted 
by  an  unknown  colored  or  non-uniform  Gaussian  noise.  A new  time  delay  estimator 
based  on  the  ML  approach  is  derived.  However,  this  estimator  requires  the  solution 
of  an  optimization  problem  whose  objective  function  is  highly  oscillatory.  Inspired 
by  related  work  [110],  a new  approach  is  proposed  to  obtain  an  initial  time  delay  es- 
timate by  assuming  that  the  amplitude  of  the  reflected  signal  is  complex-valued  and 
then  to  refine  the  initial  estimate  based  on  the  true  cost  function.  A computationally 
efficient  algorithm,  referred  to  as  the  direct  matching  method  (DMM),  is  presented 
to  obtain  the  refined  estimate.  As  shown  in  this  chapter,  our  algorithm  performs 
much  better  than  the  estimator  in  [49]  and  the  RMSEs  of  the  parameter  estimates  of 
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the  new  method  approach  the  CRBs  as  the  SNR  increases.  It  is  also  shown  that  for 
the  case  in  which  the  amplitude  is  no  longer  real-valued,  but  complex-valued  with  a 
small  phase  angle  (i.e. , in  the  case  of  a small  model  mismatch),  the  new  method  can 
still  perform  better  than  the  traditional  method  for  a low  to  moderate  SNR. 

3.2  Problem  Formulation 

Assume  that  the  transmitted  ultrasonic  signal  is  a sinusoid: 

st(t)  = sin(27r/it  + 0i),  (3.1) 

where  fx  is  the  transmitted  frequency  and  <f>i  is  the  initial  phase.  When  the  sound 
waveform  is  normally  incident  on  the  interface  between  two  media,  the  wave  is  partly 
reflected  and  partly  transmitted  into  the  second  medium.  The  reflection  coefficient 
is  defined  as  R = p~  /p+,  where  p+  is  the  pressure  of  the  incident  wave  and  p~  is  the 
pressure  of  the  reflected  wave.  This  can  also  be  expressed  as  R = (Z2  — Zi)/(Z2  + Zi), 
where  the  Zx  and  Z2  are  the  characteristic  impedances  of  the  transmission  media 
[14].  When  the  impedance  of  the  second  medium  is  much  larger  than  that  of  the  first 
medium,  for  instance,  when  the  aerial  sound  waves  impinge  onto  a water  surface,  R 
is  approximate  to  unity1,  which  implies  that  the  reflected  signal  is  an  exact  replica  of 
the  incident  signal.  In  this  case,  the  signal  received  by  the  transducer  can  be  written 
as 

8r(t)  = Pam[2irfi(t-T)  + <l>i],  (3.2) 

where  r is  the  round-trip  time  delay  between  the  sensor  and  the  reflecting  boundary 
and  ft  is  the  amplitude,  which  is  considered  to  be  a real-valued  nonnegative  unknown. 
Therefore,  the  phase-shift  between  the  transmitted  and  received  signal  due  to  the 
time  delay  is  27t/i t.  In  order  to  measure  this  phase-shift,  standard  I/Q  (Inphase  and 

Quadrature)  processing  is  performed  on  the  received  signal  [33].  This  means  that  sr(t) 

:At  STP  (1  atm,  20°C),  the  impedances  for  air  and  water  are  Z\  — 415  MKS  rayls  and  Z 2 = 
1.48  x 106  MKS  rayls,  respectively.  Thus  in  this  case,  R = 0.9994. 
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is  demodulated  with  the  inphase  and  quadrature  components  of  st(t),  sin^/if + 
and  cos(27r/it  + (j) i),  respectively.  The  demodulated  signals  are  passed  through  low- 
pass  filters  and  the  real  and  imaginary  parts  of  /dej27r^lT  are  obtained  by  measuring 
the  output  amplitudes  of  the  two  low-pass  filters. 

Assume  that  N subsequent  and  independent  measurements  are  taken  using 
the  frequency  /i.  Then  each  measurement  value  can  be  modeled  as 

x^n)  = Pej2nflT + ey{n),  n = 1, 2, ...,  N,  (3.3) 

where  we  assume  that  e\ (n)  is  a zero-mean  circular  symmetric  complex  white  (with 
respect  to  n)  Gaussian  random  variable  with  variance  af.  Assume  further  that  a 
corresponding  set  of  N independent  measurements  are  obtained  by  transmitting  a 
signal  with  the  known  frequency  f2  (without  loss  of  generality,  we  assume  /i  < /2): 

x2(n)  = + e2(n),  n = 1, 2,  • • •,  N,  (3.4) 

where  e2(n)  is  Gaussian  noise  with  variance  erf.  We  will  in  general  assume  that  the 
noise  is  non-uniform  so  that  o\^  a\.  This  modeling  accounts  for  the  possibility  that 
the  noise  level  is  different  at  different  frequencies,  a situation  that  can,  for  example, 
arise  from  hardware  artifacts.  Furthermore,  we  will  allow  for  the  possibility  that 
d (n)  and  e2(n)  are  correlated,  a situation  that  could  arise,  for  instance,  if  the  two 
measurements  (at  time  n)  at  frequency  f\  and  /2  are  taken  simultaneously  using  two 
transducers.  We  will  see  that  this  data  model  leads  not  only  to  a robust  estimation 
algorithm,  but  is  also  mathematically  efficient. 

The  data  model  can  be  expressed  using  matrix  notation  as  follows.  Define  the 
vector  of  measurements 

x(n)  = [ xi(n)  x2{n)  ]T  e C2xl,  (3.5) 

where  (-)r  denotes  the  transpose,  let 


e(n)  = [ei(n)  e2(n)  ]T  € C2*1, 


(3.6) 
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be  a vector  containing  the  noise  samples  and  let 

b(r)  = [ ei2nflT  ej2nhT  ]T. 

Then  the  data  model  in  (3.3)  and  (3.4)  can  be  written  in  a matrix  form  as 


x(n)  = /3b  + e(n),  n = 1, 2, N, 


(3.8) 


where  e(n)^=1  are  independently  and  identically  distributed  zero-mean  circularly 
symmetric  complex  Gaussian  random  vectors  with  an  unknown  positive  definite  co- 
variance  matrix 


Q 


P&1&2 

ol 


p (7\<J2  u2 

and  here  the  complex  correlation  coefficient  p is  defined  as 

E(e*e2) 


P = 


0\02 


(3.9) 


(3.10) 


Hereafter  E(x)  is  the  expected  value  of  x and  (•)*  denotes  the  complex  conjugate. 
The  problem  of  interest  in  this  chapter  is  to  estimate  {r,  /?}  from  x(n),  n = 1,2, ...,  N. 


3.3  Time  Delay  Estimation 

3.3.1  ML  Estimator 


The  log-likelihood  function  of  x(n)  is  proportional  to  (within  an  additive  con- 
stant): 

1 N 

C = — ln[det(Q)]  - tr{Q-1—  ^[x(n)  - /?b][x(n)  - (3.11) 

71=1 

where  det(-)  denotes  the  determinant  of  a matrix,  [-]ff  denotes  the  conjugate  trans- 
pose, and  tr{  • } stands  for  the  trace  of  the  matrix  [42], 

It  is  easy  to  show  that  the  matrix  Q that  maximizes  (3.11)  is  (e.g.,  [63,  25]) 

Q = ^ X^tx(n)  - ^bl[x(n)  ~ ■ (3-12) 

71  — 1 

Insertion  of  (3.12)  into  (3.11)  shows  that  maximizing  (3.11)  is  equivalent  to  minimiz- 
ing 

£i  = det  ^ J^[x(n)  - 0b][x(n)  - /3b]H  j = det(G),  (3.13) 
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where 

G = “ 5^[x(n)  - /3b][x(n)  “ Pb]H • (3-14) 

n— 1 

Let 

1 N 

x=— J^x(n),  (3.15) 

71=1 

and 

1 N 

R-xx  = x(n)x(n)H-  (3.1-6) 

71=1 

Then  G can  be  written  as 

G = Ria;  — /3bxH  — P*xbH  + \/3\2bbH 

= [fib  - x][/3b  - x]H  + Rxx  - xxff.  (3.17) 

With  this  notation,  the  cost  function  C\  in  (3.13)  can  be  written  as 

Ci  = det[RXI  — xxff  + (/3b  - x)(/3b  — x)H] 

= det(Q)det[I  + Q_1(/3b  - x)(/3b  - x)H] 

= det(Q)[l  + (/3b  - x)HQ_1(/3b  - x)],  (3.18) 

where 

Q = Rm  - xxH,  (3.19) 

and  where  we  have  used  the  fact  that  det(I  + AB)  = det(I  + BA)  whenever  the 
dimensions  of  A and  B are  conformable  [99].  Hence,  maximizing  C is  equivalent  to 
minimizing 

£2(t,  0)  = [/3b (r)  - x]HQ_1[/3b(r)  - x],  (3.20) 

which  is  a highly  nonlinear  optimization  problem. 

Before  we  proceed,  let  us  remark  on  the  following  two  facts.  First,  the  data 
model  in  (3.8),  the  related  optimization  problem  of  maximizing  (3.11)  and  its  solution 
are  similar  to  those  obtained  for  the  amplitude  and  phase  estimation  (APES)  filter 
[64],  Second,  it  is  evident  from  (3.20)  that  the  exact  ML  estimator  of  (r,/3)  reduces 
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to  a simple  NLS  cost  function  of  pre-whitened  data.  Note,  however,  that  Q in  (3.19) 
is  not  the  ML  estimate  of  Q. 

We  will  consider  three  different  cases:  /3  is  complex- valued,  ft  is  real-valued 
and  is  nonnegative,  respectively.  Consider  first  the  case  where  (3  is  an  unknown 
complex- valued  scalar.  Then  minimizing  £2(t,/3)  with  respect  to  r and  f3  yields  the 
estimates  of  f and  /3,  respectively,  as 


, , a |bHQ-1x|2 

tc  i = argmaxci(r)  = argmax  - , 

r r bwQ_1b 


A 


bHQ-1i 


Cl 


(3.21) 


(3.22) 


T=TC, 


bHQ-lb 

Next,  consider  the  case  of  a real-valued  /3.  Minimizing  £2(t,  (3)  with  respect  to  r and 
/?  yields 

, X A Re^b^Q"1*) 

Tc  2 = arg  max  c2(r)  = argmax — > .. — , (3-23) 

T T bHQ_1b 


Pc2  = 


Re(bHQ  Jx) 


(3.24) 


T=Tc2 


bHQ-1b 

where  Re(x)  is  the  real  part  of  x.  Finally,  consider  the  case  when  /?  is  a nonnegative 
unknown.  Minimizing  C2(t,/3)  for  a fixed  r with  respect  to  (3  yields 


Re(bHQ  *x)  ( Re(bHQ  xx) 

p(r)  - u 

b^Q^b  \ b^Q^b 

where  u(x)  is  the  unit  step  function: 

1 if  x > 0, 


(3.25) 


u(x) 


Insertion  of  (3.25)  into  (3.20)  gives 


fC3  — argmaxcs(r)  = argmax 


0 if  x < 0. 


Re2(bHQ  lx)  / Re(b^Q  xx) 


(3.26) 


b^Q-Jb 


-u 


bHQ-1b 


(3.27) 


Once  we  have  obtained  fC3,  the  corresponding  /3C3  is  obtained  from  (3.25)  as  /3C3  = 

m3). 
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According  to  the  parsimony  principle  [97],  if  /3  is  real- valued,  the  estimates 
obtained  by  maximizing  C2(t)  and  Cs(t)  should  be  more  accurate  than  those  obtained 
by  maximizing  c\{r)  due  to  the  different  constraints  on  j5.  Although  (3.23)  and 
(3.27)  are  simple  1-D  search  problems,  the  objective  functions  c2(r)  and  c3(r)  are 
highly  oscillatory  and  have  numerous  closely  spaced  local  maxima  which  make  it 
very  difficult  to  find  the  global  maximum.  In  the  section  below,  we  present  a fast 
algorithm  to  cope  with  these  optimization  problems. 

3.3.2  Estimation  Algorithm 


To  solve  the  optimization  problems  in  previous  section,  we  first  introduce  the 
following  notation.  Let 


7i 
v*  72 


Vv  = arg(z/), 

Xi  = — JJa^n),  * = 1,2, 

71  — 1 

i’xi  = arg(xj),  * = 1,2, 


(3.28) 

(3.29) 

(3.30) 

(3.31) 


where  Q is  defined  in  (3.19)  and  arg(x)  is  the  phase  of  x. 

Maximization  of  ci(r).  Consider  first  the  optimization  problem  for  the 
cost  function  ci(r),  which  after  some  straightforward  algebra  can  be  written  as 


A + \B\ cos  (ip  + ipB) 

7i  + 72  + 2|i/|  cos  (ip  + ipu)  ’ 


(3.32) 


where 


A = |a^i  |2  (Ti  + M2)  + l^2|2(72  + M2)  + 2Re[x1x;(7i  + 72K],  (3.33) 

B = 2[z^(|xi|27i  + |x2|272)  + Xix^i^  + x\x2v2],  (3.34) 

-0b  = arg(B),  (3.35) 


</>  = 2? r(/2  - /i)r. 


(3.36) 
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Differentiating  c^r)  with  respect  to  and  setting  the  derivative  to  zero  yields 


B | sin(^  + 7/>B)[7i  + 72  + 2M  cos (V’  + Vv)]  - [A+\B\  cos(i/)  + ^B)][2|^|  sin('0  + '0I/)]  = 0. 

(3.37) 

Let 

p = \B\(ji  +72)  cos (V>b)  - 2|i'|Acos(t/v),  (3.38) 

q = |B|(7!  + 72)sin(^B)  - 2|z/|A sin(i/v),  (3.39) 

^ = arg(p  + jg).  (3.40) 


Then  the  solution  of  (3.37)  is 


7 • / 2|i/||S|sin(V>B  lT, 

w = arcsin . ) — T. 

v/Pr+^ 

If  it  is  known  a priori  that  rmax  < l/(/2  — /1 ) , the  estimate  of  r is 


(3.41) 


'Cl 


± 

2tt(/2-/i) 

iP+2tt 

MS2-S1) 


if  %l>  > 0, 
if  -0  < 0. 


(3.42) 


Note  that  in  the  special  case  that  it  is  known  that  Q = a2I  where  I is  the  identity 
matrix  and  a2  = a2  — o\,  the  estimator  in  (3.42)  can  be  simplified  to 


V'i2_V’x1 

2tt(/2-/i) 

{)px2-)pxl)+  27T 

Mh-h) 


if  - V’l!  > 0, 
if  V>x2  - V’ai  < 0, 


(3.43) 


which  is  the  same  expression  as  for  the  traditional  phase-shift  approach  in  [49]. 

DMM.  Although  a closed-form  solution  to  maximizing  C\{r)  with  respect 
to  t can  be  found,  it  is  less  accurate  than  the  solution  obtained  by  maximizing  c2(r) 
and  c3(t),  if  p is  real- valued.  However,  it  is  difficult  to  maximize  c2(r)  and  c3(r) 
directly  due  to  the  fact  that  both  of  them  are  highly  oscillatory  with  c3(r)  being  even 
non-smooth. 

In  this  work,  we  propose  to  maximize  the  cost  function  Ci(r)  by  first  obtaining 
an  initial  estimate  fCl,  and  then  refine  it  based  on  C2(r)  or  c3(r)  respectively.  Our 
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algorithm,  which  does  not  require  a search  through  the  parameter  space  of  r,  is 
described  below  and  will  be  called  the  direct  matching  method  (DMM)  algorithm. 
Rewrite  the  cost  function  C2(r)  as 

c2(r)  =£i(r)  -gKr),  (3.44) 


where 

9i(t)  = (bHQ_1b)_1 

= [7i  +72  + 2|^|cos(27t(/2  - fi)r  + 0„)]_1,  (3.45) 

and 

g2{r)  = Re^Q-1*) 

= 71|x1|  cos(27r/ir  - 02l)  + 72 1^2 1 cos(27 r/2r  - 022) 

+ |j/||xi|  cos(27r/2r  - 02l  + Vv) 

+MN  cos(2t r/iT  - -0X2  - 0«/)-  (3-46) 

It  can  be  noted  that  the  denominator  of  gi(r)  contains  a sinusoid  with  (low)  frequency 
/2  — /1  and  that  p2(r)  is  highly  oscillatory  with  an  oscillation  frequency  somewhere 
between  f\  and  f2.  Therefore,  the  cost  function  c2(t)  is  highly  oscillatory.  Differen- 
tiating c2(r)  with  respect  to  r yields 

c2(t)  = g[{T)  ■ g2(r)  + 2gl{r)  ■ g2(r)  ■ g2{r).  (3.47) 

Next,  we  expand  c2(r)  around  fCl  to  a first-order  approximation  and  find  the  solu- 
tion to  the  equation  c'2(r)  = 0,  which  corresponds  to  the  local  maximum  of  c2(r) 
around  fCl.  A Taylor  series  expansion  of  gi(r)  around  fCl  yields  to  a second-order 
approximation 


9i(r)  &a0  + Oi(r  - fCl)  + a2(r  - fCl)2, 


(3.48) 
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where 


a i 


°0  — 5l(7")|r=fc1 

= [71  + 72  + 2M  cos(27 r(/2  - /i)fCl  + t/v)]-1, 
(r) 


(3.49) 


dr 


T=TC1 


and 


a 2 


= 4tt(/2  - /i)|i/|  sin(27r(/2  - /i)fCl  + ^„)[7i  + 72 
+2|z/|  cos(27 r(/2  - /i)fCl  + t/V)]-2, 


l5V(r) 


(3.50) 


2 dT2  r=rci 

4t t2(/2  - /1)2|^|{cos(2tt(/2  - /i)fCl  + Vv)[7i  + 72 
+2|i/|  cos(27 t(/2  - /i ) tCi  + ^)]-2  + 4|i/|  sin2(27r(/2  - /i)f, 
+Vv)[7i  + 72  + 2|i/|  cos(2t r(/2  - /i)fCl  + t/^)]“3}. 


Cl 


(3.51) 


Since  j3  is  nonnegative,  we  have  the  following  approximation  for  high  SNR: 


27T  fiTCl  « 2Zj7T  + ipSi , i = 1, 2, 


where 


2,:  = 


27T/jfCl  - ^ 
27 r 


, t = 1,2, 


(3.52) 


(3.53) 


with  L-J  denoting  rounding  to  the  nearest  integer.  Note  that  the  information  that 
P > 0 is  used  explicitly  here.  In  a small  range  around  fCl,  (3.52)  yields  the  following 
approximation 

cos[27r/i(r  - fCl)]  « cos(27r/jT  - xp£i  - 2z<7r)  « 1,  i = 1,  2,  (3.54) 

and 


sin[27r/j(r  - fCl)]  « sin(27 r/,r  - « 27r/iT  - tpSi  - 2^^,  * = 1,2.  (3.55) 
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Then  g2{r)  and  g 2(r)  can  be  expanded  in  their  first-order  approximations  around  fCl 
as 

g2(r)  ^b0  + bxr , (3.56) 

and 

^(t)  ~ d0  + dir,  (3.57) 

where 

&o  = |*i|[7i  + M «»($)  - |i/|  sin(^)(V>z2  + 2z2n)] 

+|^2|[72  + M cos(?/>)  + \u\  sin(V>)(V>21  + 2z17t)],  (3.58) 

&i  = 2tc\v\  sin('0)(/2|^i|  — fi\x2\),  (3.59) 


d0  = 27r|i/|sin(^)(/2|xi|  - /i|x2|) 

+27r/i(t/>x1  + 2z17t)[M|x2|  cos('0)  + 7i|®i|] 

+27t/2(^*2  +2z27t)[|i/||Si|cos(^)  +72|®2|]» 

(3.60) 

di  = -47r2|xi|[|^|/22cos(?/i)  +71/12]  - 47r2|x2|[M/i2cos(V0  +72/2],  (3-61) 


with 

$ = ^x2  ~ (3-62) 

Therefore,  around  fCl,  c2(r)  can  be  approximately  expressed  as 

c2(r)  « [ai  +2a2(r  - fCl)](60  + f>ir)2  + 2[a0  + ai(r  - fCl)](b0  + &iT)(d0  + dxr).  (3.63) 

Simplifying  (3.63),  neglecting  the  second-order  term  of  r — fCl,  and  setting  it  to  zero, 
we  obtain  the  refined  solution  as 


TDMM  — rcx 


dib201  + 2ao&oidoi 

2(02^0!  + + dob\doi  + a05oidi  + «i5oidoi) 


(3.64) 
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where 


and 


An 


bo  + b\fCl , 


(3.65) 


doi  — do  + difCl.  (3.66) 

As  a final  remark,  we  note  that  if  it  is  known  a priori  that  Q = a2 1,  i.e.,  71  = 72  and 
v = 0,  and  if  the  SNR  is  high,  i.e.,  \xi\  « |x2|,  (3.64)  can  be  simplified  to 

/22  ,2z2tt  + ipS2 , 


tDMM|q=(T2i 


fl  ,2Z!7r  + -0X1 

f2  + f2[  27 r/x  ’"fl  + fV  27 r/2 


) + 


L)- 


(3.67) 


Equation  (3.67)  can  be  interpreted  as  a weighted  summation  of  the  two  estimates  of 
t obtained  via  f\  and  /2  separately. 

The  steps  of  the  DMM  algorithm  are  summarized  as  follows: 


1.  Estimate  Q via  (3.19).  Invert  Q to  obtain  ji,  72  and  u (see  (3.28)). 


2.  Obtain  an  initial  estimate  of  fCl  via  (3.42). 

3.  Estimate  Z\  and  z2  via  (3.53). 

4.  Compute  a0,ai,a2,bo,bi,do,di  via  (3.49)-(3.51)  and  (3.58)-(3.62). 

5.  Obtain  the  final  estimate  tdmm  via  (3.64). 

Once  the  final  estimate  t is  obtained,  computing  the  estimate  /3  is  straightforward 
via  (3.24). 

Analysis  of  model  errors.  In  practical  measurement  scenarios,  due  to 
the  complex  measurement  environment  and  in  particular  the  non-ideal  rigid  wall 
reflection  and  inhomogeneities  of  the  transmission  media,  it  may  not  be  entirely 
correct  to  model  /3  as  a nonnegative  number.  In  fact,  it  can  be  argued  that  /3  should 
be  modeled  as  a complex- valued  number  (3  — \/3\e^0  where  ipp  is  a small  phase  angle. 
In  this  case,  the  estimate  tDmm,  which  is  obtained  by  assuming  0 is  nonnegative, 
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will  be  biased  due  to  the  model  mismatch.  It  can  be  shown  that  for  the  special  case 
Q = <t2I  in  (3.67),  the  bias  and  the  mean-squared  errors  (MSE)  of  the  estimate  tdmm 
are 


Bias(fDMM)  = E(fDMM  - t\t)  = 


fi  + fi  , 


(3.68) 


and 


MSE(tdmm)  = E[(tdmm  — t)2\t ] — 


fl  + f“2 

8»2(/i2  + /22)  ' [Mfl  + /I) 


+ 


(3.69) 


It  is  clear  from  the  above  equations  that  both  Bias(fpMM)  and  MSE(tdmm)  increase 
as  'ipp  increases. 


3.4  Numerical  Examples 


In  this  section,  we  present  several  computer  simulations  to  illustrate  the  per- 
formance of  our  proposed  algorithm.  The  performance  of  our  algorithm  is  compared 
to  the  CRB,  which  gives  the  minimum  attainable  variance  for  any  unbiased  estima- 
tor. The  CRB  for  our  problem  is  derived  in  Appendix  A.  In  the  simulations  below, 
we  use  N = 10,  /3  - 0.5,  fx  - 60  kHz,  /2  = 63  kHz,  r = 1.67  x 10-4  s (corresponding 
to  0.029  m if  the  sound  speed  is  344  m/s).  The  RMSE  of  each  estimate  is  obtained 
by  500  Monte  Carlo  trials  and  the  SNR  is  defined  as  2 /52/(cr2  + erf)- 

Figure  3.1  illustrates  the  cost  functions  Ci(r),  c2(r)  and  c3(r).  As  shown 
in  Figure  3.1(a)  and  Figure  3.1(b),  Ci(r),  c2(r)  and  c3(r)  share  the  same  global 
maximum  if  /3  is  nonnegative.  Note  that  both  c2(r)  and  c3(r)  are  highly  oscillatory, 
and  also  that  c3(r)  has  many  non-differentiable  points.  On  the  other  hand,  Ci(r) 
is  the  envelope  of  c2(r)  and  c3(r)  and  is  in  fact  quite  smooth.  It  is  expected  that 
maximizing  c2(r)  or  c3(r)  can  yield  a much  more  accurate  estimate  than  maximizing 
Ci(t),  due  to  the  sharper  peaks  of  c2(r)  or  c3(r).  Moreover,  since  c2(r)  has  more 
closely  spaced  local  maxima  than  c3(r),  it  is  intuitively  expected  that  c2(r)  is  more 
sensitive  to  noise  in  a certain  SNR  range  than  c3(r)  is.  If  ft  is  complex- valued,  as 
shown  in  Figure  3.1(c),  the  global  maxima  of  c2(r)  and  c3(r)  deviate  from  the  true 
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Figure  3.1:  Comparison  of  the  cost  functions  C\  (r),  c2(r)  and  c3(r)  for  different 
assumptions  on  the  signal  amplitude  /3.  No  noise  is  present.  True  time  delay  is  1.67  x 
10~4s.  (a)  (b)  (3  is  nonnegative,  three  cost  functions  share  same  global  maximum, 
(c)  /3  is  complex- valued,  the  global  maximum  of  c3(r)  deviates  away  from  the  true 
location  of  r. 


location  of  r (for  simplicity,  only  c3(r)  is  shown  in  this  figure).  Hence,  the  estimates 
based  on  c2(t)  or  c3(r)  will  be  biased  due  to  the  model  mismatch. 

Figure  3.2(a)  illustrates  the  RMSEs  of  the  different  methods  for  the  distance 
estimates.  The  distance  is  calculated  by  rc/2  with  c = 344  m/s.  The  noise  covariance 
matrix  used  here  is  Q = <r2I  and  five  different  methods  are  compared.  The  first  is 
called  “Cl”  and  uses  the  closed-form  solution  in  (3.42)  to  maximize  cost  function 
Ci(r).  The  second  and  the  third  methods  are  based  on  maximizing  c2(r)  and  c3(r) 
by  a 1-D  search  method,  and  are  referred  to  as  “C2”  and  “C3”,  respectively.  The 
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Figure  3.2:  Comparison  of  the  RMSEs  and  CRBs  for  estimating  (a)  distance  (rc/2 
with  c = 344  m/s),  (b)  /?  when  the  signal  amplitude  is  assumed  to  be  complex- 
valued (Cl),  real- valued  (C2)  and  nonnegative  (C3).  The  noise  covariance  matrix  is 
Q = a2l. 

1-D  search  is  performed  in  two  steps,  with  an  initial  estimate  via  (3.42),  followed  by 
a fine  search  using  the  fmin  function  in  Matlab.  The  fourth  is  “C3  (White)” , which 
is  a special  case  of  C3,  where  the  Q is  set  to  an  identity  matrix,  that  is,  we  assume 
the  noise  is  white.  The  last  method  is  the  DMM  given  by  (3.64)  ((3.67)  is  not  used 
in  any  of  our  simulations).  It  can  be  noted  that  Cl  cannot  approach  the  real  CRB 
even  when  the  SNR  goes  to  infinity.  However,  both  C2,  C3  and  DMM  approach 
the  real  CRB  when  the  SNR  increases.  Note  that  the  threshold  effect  for  C3  and 
DMM  occurs  earlier  than  for  C2,  which  is  intuitively  expected  since  both  C3  and 
DMM  use  the  information  that  (3  is  nonnegative.  Note  also  that  the  performances 
of  C3,  C3(White)  and  DMM  are  almost  the  same,  despite  the  fact  that  C3  and 
DMM  do  not  take  advantage  of  the  prior  knowledge  of  the  noise.  In  other  words,  the 
robustness  against  the  colored  or  non-uniform  noise  offered  by  our  data  model  and 
algorithm  comes  at  a negligible  cost  (or  rather  a negligible  performance  loss  in  the 
case  Q = <r2I).  The  reason  for  this  is  related  to  the  fact  that  our  estimation  problem 
is  decoupled  in  the  sense  that  the  Fisher  information  matrix  is  block  diagonal  (cf. 
Appendix  A).  Figure  3.2(b)  shows  the  corresponding  estimates  of  ft. 
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Figure  3.3:  Comparison  of  the  RMSEs  and  CRBs  for  estimating  the  distance  (tc/2 
with  c = 344  m/s)  when  the  signal  amplitude  is  assumed  to  be  complex- valued  (Cl), 
real- valued  (C2)  and  nonnegative  (C3).  The  noises  of  the  two  channels  are  correlated 
with  each  other  and  have  different  variances. 


Figure  3.3  illustrates  the  RMSEs  of  the  different  methods  for  the  distance 


0.93tfi  cr2j 


— 0.93(7!  <72j 


(Jo 


where 


estimates  when  the  noise  covariance  matrix  is  Q = 
o\lo\  = 2.  Note  that  the  performance  of  DMM  is  close  to  that  of  C3  and  approaches 
the  real  CRB  as  the  SNR  increases.  If  the  information  of  Q is  not  used,  as  shown  by 
C3  (White),  there  is  a gap  of  8 dB  from  the  real  CRB,  no  matter  how  high  the  SNR 


is. 


As  a further  illustration  of  the  robustness  of  our  algorithm,  Figure  3.4(a) 
shows  the  performance  of  the  C3  and  DMM  as  a function  of  the  absolute  value  of  the 
correlation  coefficient  p by  fixing  the  SNR=35  dB,  arg(p)  = 7r/2  and  a\  = o\.  The 
figure  shows  that  both  C3  and  DMM  are  close  to  the  real  CRB,  however,  C3  (White) 
deviates  from  the  real  CRB  when  \p\  approaches  unity  (i.e.,  when  the  noises  in  the 
two  channels  become  highly  correlated).  Figure  3.4(b)  shows  the  performance  of  the 
distance  estimators  as  a function  of  the  ratio  of  the  two  variances  o\/o\  by  fixing 
SNR=35  dB  and  p = 0.  It  can  be  noted  that  C3  and  DMM  can  still  achieve  the  real 
CRB,  whereas  C3  (White)  deviates  from  the  real  CRB  when  the  ratio  increases.  The 
interpretation  of  this  is  that  our  estimator  has  ability  to  estimate  the  noise  level  on 
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Figure  3.4:  Comparison  of  the  RMSEs  and  CRBs  for  estimating  the  distance  (rc/2 
with  c = 344  m/s)  (a)  as  a function  of  the  absolute  value  of  p with  arg(p)  = 7t/2, 
and  (b)  as  a function  of  o\)o\  with  p = 0.  The  signal  amplitude  is  assumed  to  be 
nonnegative  and  SNR=35  dB. 


the  two  different  frequencies,  and  to  combine  the  measurements  optimally  taking  the 
different  noise  levels  into  account. 

From  Figures  3.3  and  3.4,  we  can  see  that  our  new  estimator  has  the  ability  to 
adapt  to  an  unknown  noise  model.  Both  C2,  C3  and  DMM  can  achieve  the  real  CRB 
as  long  as  the  SNR  is  above  a certain  threshold.  However,  DMM  is  computationally 
more  efficient  than  C2  and  C3  since  it  does  not  require  any  search. 

Finally,  in  Figure  3.5  we  illustrate  the  performance  of  DMM  when  /3  is  complex- 
valued with  a small  phase  angle.  The  noise  covariance  matrix  is  Q = a2 1.  Due  to  the 
model  mismatch,  our  estimates  are  biased,  and  can  therefore  not  approach  the  real 
CRBs  as  the  SNR  increases.  However,  if  the  phase  angle  ipp  is  small  and  the  SNR 
is  low  to  moderate,  which  means  that  the  estimation  error  caused  by  the  noise  is 
larger  than  the  constant  bias  due  to  the  model  mismatch,  DMM  still  performs  better 
than  the  traditional  method.  When  the  SNR  increases,  the  constant  bias  caused  by 
the  model  mismatch  becomes  dominant  compared  to  the  error  caused  by  the  noise, 
and  the  traditional  method  outperforms  DMM.  Note  that  the  RMSE  of  /3  for  the 
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SNR  (in  dB)  SNR  (in  dB) 

(a)  (b) 

Figure  3.5:  Comparison  of  the  RMSEs  and  CRBs  by  using  DMM  to  estimate  (a) 
distance  (rc/2  with  c = 344  m/s)  and  (b)  /?.  Here  the  phase  angle  of  /3  is  0°,  1°,  2°, 
4°,  respectively.  The  noise  covariance  matrix  is  Q = a2 1. 

traditional  method  is  lower  than  the  corresponding  CRBs  at  the  low  SNRs,  which  is 
due  to  the  fact  that  it  is  a biased  estimator. 

3.5  Summary 

In  this  chapter,  a new  time  delay  estimator  for  the  phase-shift  approach  is 
presented.  We  consider  the  case  when  the  amplitude  of  the  signal  is  a nonnegative 
unknown  and  the  measurements  are  corrupted  by  Gaussian  noise  with  unknown  co- 
variance.  To  deal  with  the  induced  optimization  problem,  whose  cost  function  is 
highly  oscillatory,  we  propose  a method  called  DMM  that  first  finds  an  initial  esti- 
mate of  the  unknown  parameter  by  maximizing  a smoother  cost  function,  and  then 
refines  this  estimate  based  on  the  true  cost  function.  The  RMSEs  of  the  estimates 
obtained  by  using  our  new  algorithm  approach  the  corresponding  CRBs  as  the  SNR 
increases.  Since  DMM  does  not  require  any  search  over  the  parameter  space,  it  is 
suitable  for  a practical  implementation.  It  is  also  shown  that  when  the  model  mis- 
match occurs,  if  the  phase  angle  of  the  reflection  coefficient  is  small,  the  DMM  can 
still  perform  better  than  the  traditional  method  for  a low  to  moderate  SNR. 


CHAPTER  4 

MULTI-FREQUENCY  TECHNIQUE 

4.1  Introduction 


In  the  previous  chapter,  we  consider  the  phase-shift  approach  in  which  two 
frequencies  are  used.  It  is  a straightforward  step  to  consider  a more  general  case  in 
which  an  arbitrary  number  of  frequencies  can  be  used.  However,  new  problems  arise. 
For  example,  how  to  measure  the  phase-shifts  for  many  frequencies  efficiently,  and 
what  is  the  corresponding  optimal  signal  processing  approach. 

In  this  chapter,  we  consider  a measurement  scheme  based  on  a multi-frequency 
excitation  and  detection  sensor  system.  A specially  designed  excitation  waveform 
with  the  orthogonality  property  is  introduced  to  simplify  the  procedure  of  measuring 
the  phase-shift  at  each  frequency.  We  then  formulate  a general  data  model  applicable 
to  this  multi-frequency  measurement  scheme  where  the  amplitude  of  the  reflected 
signal  is  considered  to  be  a nonnegative  unknown  based  on  the  assumption  of  the 
acoustically  hard  reflection.  According  to  this  data  model,  a new  time  delay  estimator 
based  on  the  NLS  fitting  criterion  is  derived.  However,  this  estimator  requires  the 
optimization  of  a 1-D  cost  function  that  is  highly  oscillatory.  The  direct  search 
method  makes  this  estimator  computationally  too  expensive  to  be  used  for  real- 
time implementations.  By  using  an  idea  similar  to  the  one  adopted  in  the  previous 
chapter,  we  present  a two-stage  fast  optimization  algorithm  by  first  obtaining  an 
initial  estimate  via  taking  advantage  of  the  special  structure  of  the  data  model,  then 
refining  it  using  a fast-matching  method  (FMM)  based  on  the  true  cost  function. 
Simulation  results  show  that  our  method  avoids  the  search  over  the  parameter  space 
and  yet  there  are  no  obvious  performance  degradations  when  compared  to  the  direct 
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search  method.  The  numerical  results  also  show  that  the  RMSEs  of  our  new  method 
can  approach  the  corresponding  CRBs  as  the  SNR  increases. 

4.2  Data  Model  and  Problem  Formulation 


Assume  that  a transmitter  is  continuously  driven  by  K different  sinusoids,  and 
herein  we  consider  one  transmitter  working  at  K different  frequencies  simultaneously. 
In  this  case,  the  desired  excitation  signal  can  be  written  as 

K- 1 

Sd(t)  = ^2  cos(27r/jfct  + (4-1) 

fc= o 

where  and  ipk  denote  the  known  frequency  and  initial  phase  of  the  kth  sinusoid, 
respectively.  Let  fk  be  uniformly  distributed  in  the  frequency  domain  and  have  the 
form 

fk  = fo  + kAf,  A;  = 0,1 AT-  1,  (4.2) 


where  A / denotes  the  frequency  increment.  To  avoid  phase  ambiguities  in  phase 
shift  measurements  in  the  reflected  signal,  A / satisfies 

A/  < — , (4.3) 

Anax 

where  rmax  denotes  the  maximum  possible  time  delay  in  a practical  measurement 
scenario.  Let 


Tr  = 


1 

A/’ 


(4.4) 


and  we  choose  /0  to  satisfy 


foTi  = q, 


(4.5) 


with  q being  a positive  integer.  Thus  we  design  Sd(t)  so  that  it  is  a periodic  signal 
with  period  T/  and  the  K sinusoidal  signals  in  (4.1)  are  orthogonal  to  each  other 
within  Tj  as  well.  As  we  will  see,  by  taking  advantage  of  this  orthogonality  between 
the  different  frequencies,  we  can  easily  measure  the  phase  shifts  corresponding  to  all 
of  the  frequencies  at  the  same  time.  The  received  signal  can  be  expressed  as 

K- 1 

y{t)  = (/fc)l  cos{2tt fk(t  - r)  + + arg [#(/*)]}  + e(t),  (4.6) 

o 
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where  j3  and  r denote  the  unknown  amplitude  and  time  delay  respectively,  H(fk) 
represents  the  known  complex-valued  transfer  function  due  to  both  the  transmit- 
ting and  the  receiving  channels  of  the  sensor  system,  arg(rr)  is  the  phase  of  x with 
arg(x)  € [ — 7r  7 r ],  and  e(t)  denotes  the  real- valued  receiver  noise,  which  is  modeled 
as  a zero-mean  Gaussian  random  process.  (In  a practical  situation,  y(t)  may  also 
contain  a signal  leaked  directly  from  the  transmitter,  which  is  the  so-called  crosstalk. 
However,  the  crosstalk  is  independent  of  the  reflected  signal,  and  thus  it  can  be 
recorded  in  advance  and  subtracted  out  from  y(t)  thereafter.)  Note  that  the  signal 
amplitude  /3  depends  on  both  the  reflection  coefficient  and  the  transmission  medium. 
In  most  ultrasonic  distance  measurement  scenarios,  the  target  surface  is  acoustically 
much  harder  than  the  medium.  For  instance,  when  an  aerial  ultrasound  wave  im- 
pinges onto  a water  surface,  which  can  be  modeled  as  a rigid  wall  reflection  where  the 
reflection  coefficient  is  close  to  unity  [14].  Hence  /3  is  only  related  to  the  attenuation 
caused  by  the  transmission  media  and  is  considered  to  be  a nonnegative  unknown. 
In  fact,  the  attenuation  due  to  the  media  is  frequency  dependent  as  well  [14].  How- 
ever, most  ultrasonic  transducers  can  be  modeled  as  a narrowband  system,  and  thus 
for  simplicity,  we  can  assume  that  the  media  attenuation  factor  is  approximately  a 
constant  within  such  a small  frequency  range.  Then,  after  A/D  conversion,  (6)  can 
be  written  as 

K- 1 

y{nTs)  = ^ (3\H(fk)\cos{2nfk(nTs-T)+'ipk+Mg[H{fk)]}+e(nTs),n  = 0,  l,...,iV-l, 

k=0 

(4,7) 

where  N is  the  total  number  of  the  samples,  and  Ts  is  the  sampling  period  and  is 
equal  to  the  reciprocal  of  the  sampling  frequency  fs  satisfying 

/,  = 4 > 2 /«•_„  (4.8) 

s 

and 


(4.9) 
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with  M being  a positive  integer.  Note  that  (4.8)  is  required  by  the  Nyquist  sampling 
principle  and  (4.9)  is  to  guarantee  that  the  sampled  sinusoidal  signals  of  different 
frequencies  in  (4.7)  are  still  orthogonal  to  each  other  within  the  interval  of  M samples 
in  the  discrete-time  domain.  The  sampled  received  signal  can  be  further  simplified 
as 

K- 1 

y(n)  = Y,  P\HU*) I C0S{27rAn  - 2tt /*r  + *l>k  + *rg[H{fk)]}  + e(n),  n = 0, 1, ...,  N - 1, 

fc=0 

(4.10) 

where 


fk  = fkTs  = 


Q_  k_ 

T[  Tr 


Tl 

M 


q + k 
M ’ 


A;  = 0,1,...,  if-  1, 


(4.11) 


and  e(n)  is  assumed  to  be  a zero-mean  white  Gaussian  random  process  with  variance 
a2.  The  problem  of  interest  herein  is  to  estimate  the  unknown  time  delay  r from 
{y(n)}„=Q  when  the  exciting  signal  sj(t)  as  well  as  the  transfer  function  H(fk)  are 
known.  It  can  be  seen  from  (4.10)  that  the  time  delay  r is  converted  into  the  phase 
shift  of  27t fkr  for  each  sinusoid  which  can  be  measured  by  comparing  the  received 
signal  with  the  reference  signal  without  the  phase  shift.  Several  techniques  can  be 
used  to  measure  the  phase  shift,  including  the  “square  wave”  method  and  the  “zero- 
crossing” method  [61,  108,  28,  69].  However,  these  techniques  can  only  measure  the 
phase  shift  of  one  frequency  at  a time.  For  our  case,  however,  by  taking  advantage 
of  the  orthogonality  of  the  sampled  signal  for  different  frequencies  within  the  length 
M,  we  can  obtain  the  phase  shifts  corresponding  to  all  of  the  frequencies  at  the  same 
time. 


Assume  that  the  number  of  the  measurements  is  N — LM,  where  L is  a 
positive  integer.  We  divide  {^(n)}^1  into  L non-overlapping  subsequence  of  length 
M and  use  the  Ith.  subsequence,  l = 0,1, ...,  L — 1,  to  form 


y (l)  = [y(lM)  y(lM  + 1)  ...  y{lM  + M - 1)  ]r,  l = 0, 1, ...,  L - 1,  (4.12) 
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where  (-)T  denotes  the  transpose.  Define  the  vector  of  the  reference  signal  as 

a(fc)  = [ 1 ej2v^  ...  ej27r^(M_1)  ]T,  k = 0,1,  ...,K  - 1.  (4.13) 

Since  the  if  sinusoids  in  {y{n)}^Q  are  orthogonal  to  each  other  within  each  subse- 
quence, (see  (4.4),  (4.5)  and  (4.9)),  the  inner  product  between  y (/)  and  a(k)  has  the 
form 

xt(k)  = aH{k)y(l)  = pW{k)e~j2*fkT  + ut(k),  k = 0, 1, ...,  K - 1,  (4.14) 

where  (-)H  denotes  the  conjugate  transpose, 

W(k)  = y/f(/*)^‘,  A = 0,1, ...,  if  - 1,  (4.15) 

and  ui(k ) is  the  noise  term.  Note  that  the  above  computation  can  be  regarded  as 
obtaining  the  discrete-time  Fourier  transform  (DTFT)  of  y ; at  frequencies  /0  to  /k-i, 
which  can  be  easily  implemented  by  currently  available  dedicated  high  performance 
DSP  chips.  Since  {a(A:)}j^_01  are  orthogonal  to  each  other,  {iq (&)}&= q1  *s  zero' 


mean  circular  symmetric  white  Gaussian  noise.  Define 

x(()  = [*,(0)  i,(l)  ...  x,(K  - 1)  ]T  6 CK*',  (4.16) 

W = diag{W(0),  W(l) W(K  - 1)},  (4.17) 

b(r)  = [ e~j27cfoT  e~j2nflT  ...  e~j27TfK-lT  }T  (4.18) 

_ [ g-j27r/0r  e-j27r(/0+A/)r  . , . e-j2n[f0+(K-l)Af]r  j T ^ 

u(()  = [n,(0)  u,(l)  ...  u,(K  - 1)  ]T  € CK*\  (4.19) 

and 

d(r)  = Wb(r).  (4.20) 

Then,  (4.14)  can  be  written  in  a more  compact  form  as 

x(O  = 0d(T)  + u(O,  / = 0, 1, ...,  L — 1.  (4.21) 


The  goal  here  is  to  estimate  the  unknown  time  delay  r from  the  data  set  {x(/)}^=01. 
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4.3  Parameter  Estimation 

4.3.1  NLS  Estimator 

We  consider  below  estimating  the  unknown  parameters  by  minimizing  the 
following  NLS  criterion 

A(r,/5)  = xil|x(()-/3d(r)||2,  (4.22) 

1=0 

where  ||  • ||  denotes  the  Euclidean  norm.  Note  that  for  the  white  noise  case,  the  above 
NLS  approach  is  the  same  as  the  ML  method. 

It  can  be  shown  that  the  cost  function  (4.22)  is  equivalent  to 

£2(r,«=||x-/?d(r)||2,  (4.23) 

where 

1 L_1 

*=t£x(  l).  (4.24) 

^ 1=0 

We  consider  next  the  nonlinear  optimization  problem  of  minimizing  (4.23). 

Assume  that  /3  is  nonnegative,  then  the  cost  function  £2  can  be  expressed  as 

C2  = [x  - /9d(r)]H[x  - /3d(r)]  (4.25) 

= xHx  — /3xHd(r)  — /ddH(r)x  + id2di/(r)d(r) 

= xHx  4-  /32dH  (T)d(r)  - 2/3Re[dH(r)x], 

where  Re(a:)  is  the  real  part  of  x.  Since  x^x  + /?2dH(r)d(r)  is  independent  of  r 
and  /?  is  nonnegative,  the  r that  minimizes  £2  is  the  one  that  maximizes  the  term 
2/3Re[d/7(r)x],  which  has  the  form 

fCl  = argmaxci(r)  = argmaxRe[d//(r)x].  (4-26) 

t r 

Note  that  the  cost  function  ci(r)  can  be  expressed  as 

Ci(t)  = Re[d//(r)x] 

= Re[bH(r)(WHx)] 

K- 1 

= ^2  |ffe|  cos(27 xfkT  - 6k), 

fe= 0 


(4.27) 
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where 

xk  = W*{k)x(k),  k = 0,1,...,  K - 1,  (4.28) 

and 

8k  = -axg(xk),  k = 0,  1.  (4.29) 

In  the  absence  of  noise, 

0k  = 27r/fcT  - 22fc7r,  A:  = 0, 1, ...,  if  — 1,  (4.30) 

with  being  integers.  It  can  be  seen  from  (4.27)  that  Ci(r)  is  a sum  of  M 

cosine  functions  of  r with  the  frequencies  /0  to  Jk-x-  Therefore,  Ci(r)  is  highly 
oscillatory,  which  makes  it  very  difficult  to  find  its  global  maximum  directly. 

If  /?  were  complex-valued,  it  can  be  shown  that  minimizing  £2  with  respect  to 
r would  yield 

fc,  = argmaxc2(r)  = argmax  |dH(r)x|  . (4-31) 

T T 1 

According  to  the  parsimony  principle,  if  /3  is  indeed  nonnegative,  the  estimate  ob- 
tained by  maximizing  Ci(r)  should  be  more  accurate  than  that  obtained  by  maxi- 
mizing C2(r).  Figure  4.1  shows  the  two  normalized  cost  functions  in  the  absence  of 
noise.  It  can  be  seen  that  the  true  cost  function  Ci(r)  denoted  as  Cl  with  solid  line  is 
highly  oscillatory,  and  therefore,  is  very  difficult  to  be  optimized  directly.  However, 
C2(r),  which  is  obtained  by  assuming  (3  being  complex- valued  and  is  denoted  as  C2 
with  dotted  line  in  the  figure,  is  the  envelope  of  Ci(r)  and  is  quite  smooth.  Note 
also  that  the  two  cost  functions  share  the  same  global  maximum  at  the  true  value  of 
r due  to  the  absence  of  noise.  It  is  evident  from  this  figure  that  maximizing  Ci(r) 
can  yield  a much  more  accurate  estimate  of  r due  to  its  sharper  dominant  peak  than 
maximizing  C2(r).  Below,  we  present  a two-stage  fast  algorithm  to  optimize  the  true 
cost  function  Ci(r). 
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Figure  4.1:  Comparison  of  the  cost  functions  Ci(t)  and  c2(r)  for  different  assumptions 
on  the  signal  amplitude  /3.  The  true  time  delay  is  0.1647  ms.  Cl  (solid  line)  denotes 
the  true  cost  function  Ci(r),  which  is  obtained  by  assuming  /3  being  nonnegative. 
C2  (dotted  line)  denotes  the  cost  function  c2(r),  which  is  obtained  by  assuming 
being  complex-valued.  The  dash-dot  line  indicates  the  true  value  of  r.  The  two  cost 
functions  share  the  same  global  maximum  in  the  absence  of  noise. 


4.3.2  Initial  Estimate  of  r 


To  optimize  the  cost  function  c^r),  which  is  highly  oscillatory,  it  is  reasonable 
for  us  to  obtain  an  initial  estimate  of  r followed  with  a refinement  step.  One  possi- 
ble initial  estimation  method  is  to  optimize  the  cost  function  c2(r),  which  is  much 
smoother  compared  with  Ci(r).  If  K = 2,  it  can  be  shown  that  there  exists  a closed 
form  solution  for  maximizing  c2(r).  However,  for  K > 2,  there  is  no  closed  form  solu- 
tion and  we  could  resort  to  a search-based  optimization  approach.  We  could  use  FFT 
with  zero-padding  to  locate  a dominant  peak  and  then  perform  a fine  search  nearby 
the  peak  location  [110].  Instead  of  maximizing  c2(r)  based  on  a search  method,  we 
propose  to  obtain  an  initial  estimate  of  r via  a closed  form  solution  by  exploiting  the 
special  structure  of  the  data  model  in  (4.14)  by  optimizing  a different  cost  function 
(see  (4.41)  later  on).  We  will  show  that  this  method  gives  comparable  performance 
to  the  search  method  but  is  computationally  much  simpler. 

Note  that  in  the  absence  of  noise,  the  data  model  in  (4.14)  can  be  written  as 


xi(k)  = e jAe 


W(k) 
W(k  - 1) 


xt(k  - 1),  k - 1,2,  ...,K  - 1, 


(4.32) 
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where 


A d = 27tA/t. 


In  the  presence  of  noise,  we  have 

_ -A.  W(k) 


xi(k)  = e 


W(k  - 1) 


xt(k  - 1)  + ui{k),  k = 1,2,  ...,K  - 1, 


where 


ut(k)  = ui(k)  - e jAe  W{k) 


-ui(k  - 1). 


(4.33) 


(4.34) 


(4.35) 


W{k  - 1) 

The  problem  of  interest  herein  is  to  estimate  Ad  based  on  the  measurements  xi(k). 
Define 

Zi{l)  = [xi(l)  xt{ 2)  ...  xi(K  - 1)  ]T,  (4.36) 

Xo(l)  = [xi(0)  xi(  1)  ...  xt(K  - 2)  ]r,  (4.37) 

u(/)  = [«/(!)  txi(2)  ...  u/(^-l)]T,  (4.38) 


and 


W = diag 


W{  1)  W(2)  W(K-l) 


(4.39) 


IT(0)’  W( l)’-"’  W{K-  2)  j ‘ 

Then  (4.34)  can  be  written  in  a matrix  form  as 

xx(/)  = e‘jAJWx0(l)  + u(/),  l = 0,1,...,  L-  1.  (4.40) 

The  unknown  Ad  can  be  estimated  by  minimizing  the  following  NLS  criterion 

L- 1 


£s(A0)=£  MO-  e-*“W*o(l) 


(4.41) 


/=0 


After  some  straightforward  algebra,  the  estimate  Ad  is  obtained  as  follows 


Ad  = — arg 


L- 1 


^x0»(l)W»x1(0 


z=o 


(4.42) 


Once  Ad  is  available,  the  initial  estimate  of  r is  obtained  by 

if  ^°’ 

A0+2t r 
2ttA  / 


TO  = 


if  Ad  < 0. 


(4.43) 
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4.3.3  Refined  Estimate  of  r 

We  now  refine  f0  based  on  the  true  cost  function  ci(r).  We  present  below  a 
fast  optimization  method,  referred  to  as  the  fast-matching  method  (FMM),  to  refine 
to- 

Differentiating  ci(r)  with  respect  to  r yields 

K- 1 

ci(T)  = - l^fc|27r/fc  sin(27r/fc7-  - 9k).  (4.44) 

k=0 

The  goal  here  is  to  expand  (4.44)  around  f0  with  first-order  approximation  and  find 
the  solution  of  c\(t)  — 0,  which  corresponds  to  the  local  maximum  of  Ci(r)  around 
f0.  From  (4.30),  we  have  the  following  approximation  at  high  SNR 

Znfkto  ~ 2zkTc  + dk,  k = 0, — 1,  (4.45) 


where 


zk 


fkr0  - 9k 
27T 


, fc  = 0, 1, ...,  K — 1, 


(4.46) 


with  L-J  denoting  rounding  to  the  nearest  integer.  In  a small  range  around  f0,  using 
(4.45)  yields  the  following  approximation 

sin(27 xfkT  - 4)  = sin(27r fkr  - 9k  - 2zkn)  « sin[27r/fc(r  - f0)]  « 2nfkT  -9k-  2zk'K. 

(4.47) 

Therefore  c)  (r)  can  be  approximately  expressed  around  f0  as 

K- 1 


Ci(t)  w - 2n\xk\fk{2TrfkT  - 9k  - 2zkn). 


(4.48) 


fc=0 


Simplifying  (4.48)  and  setting  it  to  zero,  we  obtain  the  refined  solution  of  r as 


EkJ  \%k\fk{9k  + 2zkn) 


T = 


Ek=o^\f2k 


(4.49) 


This  estimator  can  be  interpreted  as  a weighted  sum  of  K estimates  of  r obtained 
from  frequencies  /o  to  fx-\  separately.  Note  that  the  weight  is  proportional  to  fk, 
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which  means  that  the  estimates  obtained  at  the  higher  frequencies  are  weighted  more 
than  those  obtained  at  the  lower  frequencies.  This  interpretation  is  intuitive  because 
the  higher  the  frequency,  the  shorter  the  wavelength,  and  hence  the  better  the  time 
delay  measurement  accuracy. 

The  completed  steps  of  our  two-stage  fast  estimation  algorithm  are  summa- 
rized as  follows: 

1.  Obtain  the  initial  estimate  of  f0  via  (4.42)  and  (4.43). 

2.  Obtain  xk,  9k  and  zk,  k = 0, 1, ...,  K - 1,  via  (4.28),  (4.29),  (4.46),  respectively. 

3.  Obtain  the  refined  t by  (4.49). 

4.4  Numerical  Examples 

We  present  below  several  numerical  examples  to  illustrate  the  performance  of 
our  proposed  algorithm.  The  performance  of  our  algorithm  is  compared  to  the  CRB, 
which  is  derived  in  Appendix  B. 

In  the  simulations  below,  we  use  K = 8,  /0  = 60  kHz,  A / = 1.5  kHz.  The 
period  of  the  signal  is  Tj  = 0.6667  ms  and  the  sampling  frequency  is  fs  = 192  kHz. 
Thus  the  number  of  the  sampled  data  per  period  is  M — 128.  Totally  L — 4 periods 
are  sampled  and  the  total  number  of  data  samples  is  N = LM  — 512.  The  initial 
phases  for  each  sinusoid  are  { 0,  0,  0,  7r,  0,  n,  7r,  0 }.  The  H(fk)  is  simulated 
based  on  lumped  modeling  as  [96,  36]:  H(fk ) = jjkB*-fj+p  > ^ = 0,1, ...,  K — 1,  with 
Bw  = 10.5  kHz  being  the  3 dB  bandwidth  of  the  transfer  function  and  fc  = 65.25 
kHz  being  the  resonant  frequency.  The  true  time  delay  is  r = 0.1647  ms,  which 
corresponds  to  0.028  m for  the  sound  speed  c = 344  m/s.  The  SNR  is  defined  as 
K/32/a2,  and  the  RMSE  of  each  estimate  is  obtained  through  200  Monte-Carlo  trials. 

Figure  4.2  shows  the  waveform  of  the  signal  applied  to  the  transmitter  and 
the  received  signal  by  the  receiver  in  the  time  domain  in  the  absence  of  noise.  Two 
periods  are  illustrated.  Note  that  by  properly  choosing  the  initial  phase  for  each 


52 


2 
1.5 
1 

o> 

-O 

a 
E 
< 


” 0 0.2  0.4  0.6  0.8  1 1.2 

Time  (Second)  x 10'3  Time  (Second)  x 10‘3 

(a)  (b) 

Figure  4.2:  Input  and  output  signals  of  the  CW-based  sensor  system  in  the  absence 
of  noise.  Two  periods  are  illustrated,  (a)  The  signal  applied  to  excite  the  transducer, 
(b)  Received  signal  with  time  delay  r = 0.1647  ms. 

exciting  sinusoid,  the  signal  energy  is  evenly  distributed  in  the  time  domain  and  the 
peak-to-average  ratio  of  the  transmitting  signal  is  efficiently  reduced,  which  decreases 
the  burden  for  both  the  transmitter  and  the  receiver  at  any  time  instant.  Note  also 
that  the  waveform  shape  of  the  received  signal  is  somewhat  changed  due  to  the  effect 
of  the  transfer  function  H(fk). 

Figure  4.3  shows  the  performance  of  our  two-stage  fast  time  delay  estimation 
algorithm.  The  distance  is  calculated  by  rc/2  with  c = 344  m/s.  First,  the  perfor- 
mances of  two  different  initial  estimation  methods  are  compared.  One  is  denoted  as 
C2(s),  which  is  based  on  optimizing  the  cost  function  c2(r)  by  the  search  method. 
To  optimize  C2(r),  we  use  FFT  with  zero-padding  (to  1024  points)  to  find  the  loca- 
tion of  the  dominant  peak,  and  then  perform  a fine  search  nearby  the  peak  location 
by  using  the  function  fmin  in  Matlab.  The  other  initial  estimation  method  is  our 
initial  fast  method,  referred  to  as  IFM.  It  can  be  seen  that  the  performance  of  C2(s) 
is  slightly  better  than  that  of  IFM  and  the  difference  between  these  two  methods 
decreases  as  the  SNR  increases.  Note  also  that  neither  of  them  can  approach  the 
true  CRB  no  matter  how  high  the  SNR  is.  Next,  according  to  the  initial  estimates, 
the  refinement  step  is  performed  based  on  the  true  cost  function  Ci(r).  Two  different 
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Figure  4.3:  RMSEs  of  the  distance  estimate  (rc/2  with  c = 344  m/s)  for  the  different 
algorithms  compared  with  the  CRB.  IFM  denotes  the  initial  estimate  obtained  by 
using  our  fast  method.  C2(s)  denotes  the  initial  estimation  by  maximizing  c2(r) 
with  the  search  method.  FMM  denotes  the  estimate  by  using  our  fast-matching 
method  based  on  the  initial  estimate  obtained  by  IFM.  Cl(s)  denotes  the  estimate 
by  maximizing  Ci(r)  with  the  search  method  based  on  the  initial  estimate  obtained 
by  C2(s). 

methods  are  compared  here  as  well.  One  is  denoted  as  Cl(s),  which  uses  fmin  to 
perform  a direct  search  on  Ci(r)  nearby  the  initial  estimate  obtained  by  C2(s).  The 
other  is  our  FMM  method,  which  refines  the  initial  estimate  obtained  by  IFM.  Note 
that  Cl(s)  shows  a slightly  lower  SNR  threshold  (about  2 dB)  than  FMM.  However, 
as  the  SNR  increases,  there  is  no  obvious  difference  between  these  two  methods  and 
both  of  them  achieve  the  true  CRB. 

Finally,  Figure  4.4  shows  the  CRBs  of  the  distance  as  well  as  the  RMSEs  of 
our  new  two-stage  fast  algorithm  as  a function  of  the  number  of  exciting  frequencies 
K when  the  SNR  is  fixed  at  6 dB.  Three  different  initial  frequencies  (i.e.,  (a)  /0  = 30 
kHz,  (b)  /o  = 45  kHz,  and  (c)  /0  = 60  kHz)  are  considered,  which  represent  the 
three  different  working  frequency  bands.  The  transfer  function  H(fk)  is  set  to  1 
for  easy  comparison.  All  of  the  other  simulation  parameters  are  the  same  as  those 
used  in  Figure  4.3.  It  can  be  seen  that  the  CRBs  decease  slightly  as  K increases, 
which  is  evident  because  the  greater  the  number  of  frequencies,  the  larger  the  signal 
bandwidth  becomes,  and  hence  the  better  the  estimation  accuracy.  It  is  worth  noting 
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the  threshold  effect  of  K where  the  RMSEs  of  the  estimates  deviate  from  the  CRBs 
as  K decreases.  This  is  due  to  the  fact  that  the  smaller  K becomes,  the  wider  the 
peak  of  the  envelope  of  the  true  cost  function,  and  hence  the  harder  it  is  to  resolve 
the  ambiguity  of  the  true  cost  function  [40].  Note  also  that  the  higher  the  working 
frequency  band  chosen,  for  instance,  with  /0  = 60  kHz,  the  lower  the  CRB  can  be 
obtained.  This  is  because  as  the  wavelength  of  the  transmitted  signal  decreases,  the 
measurement  accuracy  increases.  However,  given  the  same  number  of  frequencies 
and  the  frequency  increment,  the  higher  the  working  frequency  band,  the  higher  the 
SNR  is  needed  to  achieve  the  corresponding  CRB.  This  is  intuitive  because  by  given 
the  signal  bandwidth,  the  higher  the  center  frequency,  the  more  oscillatory  the  true 
cost  function  is  (see  (4.27)),  and  hence  the  higher  SNR  value  needed  to  resolve  the 
ambiguity  of  the  true  cost  function  [113]. 

4.5  Summary 

In  this  chapter,  a multi-frequency  CW-based  ranging  technique  is  presented 
for  proximity  ultrasonic  distance  sensors.  We  have  considered  the  case  where  the 
amplitude  of  the  received  signal  is  nonnegative  due  to  acoustically  hard  reflection. 
A two-stage  fast  optimization  algorithm  is  proposed  to  determine  the  distance  effi- 
ciently. The  RMSEs  of  the  estimates  obtained  by  using  our  method  can  approach 
the  corresponding  CRBs  as  the  SNR  increases.  Since  the  new  method  does  not  re- 
quire any  search  over  the  parameter  space,  it  is  very  suitable  for  practical  real-time 
implementations. 
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Figure  4.4:  CRBs  of  the  distance  (re/ 2 with  c = 344  m/s)  and  RMSEs  of  our  two- 
stage  estimation  algorithm  (IFM+FMM)  as  a function  of  the  number  of  frequencies 
K when  SNR=6  dB  and  H(fk ) = 1,  k = 0,1,...,  if  - 1.  Three  different  initial 
frequencies  are  considered:  (a)  /o  = 30  kHz,  (b)  /o  = 45  kHz,  and  (c)  /o  = 60  kHz. 


CHAPTER  5 
PEARS  SCHEME 

5.1  Introduction 

In  previous  chapters,  two  CW  proximity  ranging  techniques  (i.e.,  the  phase- 
shift  method  and  the  multi-frequency  technique)  have  been  proposed.  The  key  step 
for  CW-based  methods  is  to  measure  the  phase-shift  for  each  transmitted  frequency. 
In  Chapter  4,  to  simplify  the  phase-shift  measurement  procedures,  the  specially  de- 
signed CW  waveform  is  presented.  At  this  point,  we  note  that  the  multi-frequency 
technique  can  be  further  extended  to  a more  general  case  where  the  transmitted  wave 
is  only  required  to  be  periodic.  This  is  because  that  after  applying  FFT  to  one  pe- 
riod of  the  sampled  signal  with  an  arbitrary  waveform,  in  the  frequency  domain,  a 
data  model  similar  to  that  used  by  the  multi-frequency  technique  can  be  formulated. 
Hence,  the  two-stage  fast  time  delay  estimation  algorithm  developed  for  the  multi- 
frequency technique  can  be  applied  to  this  new  case  with  only  minor  modifications. 
The  ability  to  deal  with  the  time  delay  estimation  problem  for  an  arbitrary  trans- 
mitted waveform  provides  more  flexibility  in  designing  a practical  proximity  ranging 
system. 

In  this  chapter,  we  first  establish  a general  data  model  for  arbitrary  trans- 
mitted waveforms.  Using  a methodology  similar  to  the  one  adopted  in  the  previous 
cases,  a new  time  delay  estimator  based  on  the  NLS  fitting  criterion  is  derived.  To 
deal  with  the  optimization  of  a 1-D  cost  function  that  is  highly  oscillatory,  a two- 
stage  fast  time  delay  estimation  algorithm,  referred  to  as  the  PEARS  (Parameter 
Estimation  for  Acoustic  Ranging  Systems)  algorithm,  is  presented.  Next,  a numeri- 
cal example  is  used  to  illustrate  the  performance  of  PEARS.  Numerical  results  show 
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that  the  RMSEs  of  the  estimates  obtained  by  using  PEARS  can  achieve  the  CRB  as 
the  SNR  increases.  Finally,  the  efficiency  of  PEARS  when  used  with  the  specially 
designed  CW  waveform  is  demonstrated  by  performing  experiments  using  commer- 
cially available  transducers.  Experimental  results  show  that  PEARS  can  be  used  to 
obtain  a high  ranging  accuracy  and  a parameter  update  rate  of  around  1.5  kHz  for  a 
distance  of  less  than  100  mm. 

5.2  Data  Model  and  Problem  Formulation 

Since  we  must  continuously  monitor  the  boundary  changes,  we  use  periodic 
waveforms  and  the  parameter  estimates  are  updated  from  one  period  to  another. 
Assume  that  a periodic  signal  with  period  T is  used  to  continuously  excite  the  trans- 
mitter as  shown  in  Figure  5.1(a).  Then  the  transmitted  signal  received  by  the  receiver 
directly  is  also  a periodic  signal  with  period  T as  shown  in  Figure  5.1(b).  Note  that 
the  waveform  in  Figure  5.1(b)  is  different  from  the  excitation  waveform  in  Figure 
5.1(a)  due  to  the  impacts  of  the  frequency  responses  of  both  the  transmitter  and  the 
receiver.  The  reflected  signal  from  a planar  boundary  is  a replica  of  the  transmitted 
signal,  and  is  shown  in  Figure  5.1(c).  To  avoid  ambiguity,  the  period  T should  satisfy 
T > where  r denotes  the  distance  between  the  centers  of  the  trans- 

mitter and  receiver,  dmax  represents  the  maximum  distance  between  the  hull  surface 
to  the  air- water  boundary  (see  Figure  1.1(b)),  and  c is  the  speed  of  sound. 

Then  the  measured  data  for  one  period  has  the  form: 

y(t)  = fis(t  — t)  + e(t),  0 <t<T,  (5.1) 

where  s(t),  0 < t < T,  represents  the  real- valued  transmitted  signal  and  T is  the 
signal  duration;  y(t ) denotes  the  real-valued  received  signal;  /3  denotes  the  amplitude 
of  the  received  signal  and  is  assumed  to  be  nonnegative  due  to  the  acoustically  hard 
boundary;  r denotes  the  time  delay;  e(f)  is  the  real-valued  noise,  which  is  modeled 
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Figure  5.1:  Diagram  of  the  signals  used  in  the  acoustic  ranging  system,  (a)  Periodic 
excitation  signal  for  the  transmitter,  (b)  Transmitted  signal  received  by  the  receiver 
directly,  (c)  Reflected  signal  from  a planar  boundary  with  the  time  delay  t and  the 
amplitude  /3  relevant  to  the  signal  in  (b). 

as  a zero-mean  Gaussian  random  process.  Note  that  s(t)  is  an  arbitrary  but  known 
waveform. 

The  received  signal  after  sampling  and  A/D  conversion  has  the  form: 

y(nTs)  = ps(nTs  - r)  + e(nTs),  n = 0, 1, N - 1,  (5.2) 

where  Ts  denotes  the  sampling  period  and  is  equal  to  the  reciprocal  of  the  sampling 
frequency  fs.  The  transmitted  signal  s(nTs)  is  known  a priori  and  {/?,r}  are  un- 
known. The  problem  of  interest  here  is  to  estimate  the  unknown  time  delay  r from 
{y(nTs)}Z~0\ 

Let  Y(m),  S(m),  and  E(m ),  m = —N/2,  —N/2  + 1, ...,  N/2  — 1,  denote  the 
discrete  Fourier  transforms  (DFTs)  of  y(nTs),  s(nTs),  and  e(nTs),  respectively.  Pro- 
vided that  the  aliasing  due  to  sampling  is  negligible,  Y (m)  can  be  written  as 

• 27 TT 

Y(m)  = -f-  E(m).  (5.3) 

We  prefer  to  work  on  the  data  model  in  the  frequency  domain  since  we  can  easily 
deal  with  the  fractional  time  delays.  Moreover,  the  resulting  parameter  estimation 
algorithm  can  be  performed  via  the  efficient  FFT  operations. 
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We  consider  below  estimating  the  unknown  parameters  by  minimizing  the 
following  NLS  criterion: 


7V/2 — 1 

£i({/3,t})=  \y (m)  - ps(m)e-’*%m 


(5.4) 


m=-N/2 

Since  both  the  transmitted  signal  s(t)  and  the  received  signal  y(t ) are  real-valued, 
their  Fourier  transforms  are  conjugate  symmetric.  Hence  the  above  cost  function  is 
equivalent  to 


A>({/3,t})=  ^2  D2(m)  Y(m)  - pS(m)e 

m=-N/  2 


(5.5) 


where  (£>(m)  = 1 }mL_w/2+1  and  D(-N/2)  = D( 0)  = j-. 
Let 


x = [ D(-N/2)Y(-N/2)  D(—N/2  + l)Y(—N/2  + 1)  •••  £>(0)F(0)  ]T,  (5.6) 

where  K = N/2  + 1 is  the  dimension  of  the  vector  x and  (-)T  denotes  the  transpose. 
Let  W be  a K x K diagonal  matrix: 

W = d\ag{D(—N/2)S(—N/2),  D{-N/2  + l)S(-N/2  + !),-••,  D(0)S(0)}.  (5.7) 


Let 

b(r)  = [e-j^~N/2)  ...  1]T 

— [ g-j27r/0T  g— j27r/ir  . . . e-j2nfK-iT  j T ^ 

where 

fk  = — 1-  + -^-,  A:  = 0, 1, ...,  K — l. 

J 2 Ts  NTS' 

Let 

d(r)  = Wb(r). 

Then  (5.5)  can  be  expressed  as 


(5.8) 


(5.9) 


(5.10) 


Aft/i,  t})  = ||x  -£d(r)||2, 


(5.11) 
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where  ||  • ||  denotes  the  Euclidean  norm.  Note  that  (5.11)  is  similar  to  the  cost  function 
in  (4.22)  which  is  derived  for  the  multi-frequency  case. 

For  the  white  noise  case,  the  above  NLS  approach  is  the  same  as  the  ML 
method.  When  the  noise  is  colored,  the  NLS  approach  can  still  give  very  accurate 
parameter  estimates  [65].  Note  that  minimizing  (5.11)  is  a highly  nonlinear  optimiza- 
tion problem. 

We  remark  that  the  data  model  in  (5.3)  and  the  cost  function  in  (5.11)  are 
similar  to  those  considered  in  [110].  However,  the  work  in  [110]  did  not  consider 
the  constraint  that  /3  is  nonnegative  and  the  proposed  optimization  algorithm  in 
[110]  involves  the  search  over  the  parameter  space.  We  present  below  an  efficient 
optimization  method,  which  is  referred  to  as  the  PEARS  (Parameter  Estimation  for 
Acoustic  Ranging  Systems)  algorithm,  to  minimize  (5.11)  by  taking  into  account  that 
/3  is  nonnegative. 


Based  on  the  derivation  similar  to  that  in  Section  4.3,  it  can  be  shown  that 
the  r minimizing  £2  with  the  nonnegative  assumption  on  f3  has  the  form: 


5.3  The  PEARS  Algorithm 


argmaxcRr)  = argmaxRe[dH(r)x[. 

T T 


(5.12) 


Note  that  the  cost  function  Ci(r)  can  be  expressed  as 


Ci(t)  = Re[dH(r)x] 


(5.13) 


K- 1 


= |£fc|cOs(27T/fcT  - ek), 


k= 0 


where 


Xk  = W*(k)x(k),  k = 0, 1, ...,  K - 1, 


(5.14) 


with  (•)*  being  the  complex  conjugate,  and 


0k  = ~ arg(x*),  k = 0, 1,...,  K - 1, 


(5.15) 
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with  arg(x)  being  the  phase  of  x.  In  the  absence  of  noise, 

9k  = 2-k fkr  - 2irzk,  k = 0, 1, K - 1,  (5.16) 

with  {zk}k=Q  being  some  integers.  It  can  be  seen  from  (5.13)  that  C\(t)  is  a sum  of 
K cosine  functions  of  r with  frequencies  /0  to  Jk-u  and  hence  is  highly  oscillatory, 
which  makes  it  very  difficult  to  find  the  global  maximum  directly. 

If  /3  were  complex- valued,  it  can  be  shown  that  minimizing  £2  with  respect  to 
t would  yield 

fC2  = argmax c2(r)  = argmax  |dH(r)x| . (5-17) 

r r 

Note  that  c2(r)  is  much  smoother  than  cost  function  C\(t).  In  Chapter  4,  we  proposed 
to  use  the  initial  fast  method  (IFM)  to  obtain  the  initial  estimate.  IFM  has  the  closed 
form  solution  for  the  initial  estimate.  However,  IFM  requires  the  division  operation 
between  two  adjunctive  spectral  measurements.  For  the  arbitrary  waveform  case,  we 
should  not  do  so  because  W(k)  could  be  zero  for  some  k.  Here  we  directly  focus 
on  maximizing  cost  function  c2(r).  Note  that  c2(r)  can  be  efficiently  maximized  by 
applying  FFT  with  zero-padding  to  WHx.  (In  our  examples,  we  zero-pad  WHx  to 
a data  length  of  4 N to  find  the  initial  estimate,  which  is  sufficiently  accurate  for  the 
second  estimation  stage.)  Once  the  initial  estimate,  say  fo,  is  available,  we  can  refine 
it  based  on  the  true  cost  function  Ci(r),  in  which  the  second  stage  of  the  algorithm 
for  the  multi-frequency  technique  can  be  employed. 

The  steps  of  the  PEARS  algorithm  can  be  summarized  as 

Stage  1: 

Obtain  the  initial  estimate  f0  of  r via  maximizing  c2(r)  by  using  zero-padding 
FFT. 

Stage  2: 

(a)  Obtain  xk  and  9k,  k = 0, 1, ...,  K — 1,  via  (5.15)  and  (5.16),  respectively. 
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(b)  Obtain  Zk  as  Zk  = 
to  the  nearest  integer. 


2nfkT0-6k 

2?r 


k — 0,1, K — 1,  with  [-J  denoting  rounding 


(c)  Obtain  the  refined  estimate  t of  r as  t = jffj { k . 

We  remark  that  the  PEARS  algorithm  is  computationally  quite  efficient.  The 
major  computational  cost  occurs  at  the  first  stage  due  to  calculating  the  initial  esti- 
mate of  fo  by  using  zero-padding  FFT.  However,  FFT  can  be  efficiently  implemented 
by  using  current  highly  dedicated  DSP  chips,  such  as  ADSP-21160  and  TMS320C6x. 
Moreover,  to  continuously  monitor  a practical  environment,  if  there  is  only  a very 
small  distance  change  between  two  consecutive  observations,  we  need  to  use  both 
stages  of  PEARS  for  only  the  first  few  parameter  estimations.  Thereafter,  the  current 
refined  estimate  f can  be  directly  used  as  the  initial  estimate  for  the  next  estimation, 
and  thus  the  first  stage  of  PEARS  can  be  totally  omitted.  Hence,  for  this  case,  the 
major  computational  cost  only  depends  on  the  second  stage  of  PEARS. 


5.4  Numerical  Example 


To  demonstrate  the  estimation  accuracy  of  the  proposed  PEARS  algorithm, 
we  present  a numerical  example  below  based  on  a windowed  chirp  signal.  The  trans- 
mitted signal  for  one  period  is  defined  as 


s(t ) = w(t)  cos 


27T  fct  + 7 


0 < t < T, 


(5.18) 


where  fc  denotes  the  carrier  frequency,  7 represents  the  chirp  rate,  and 

( 0.5  — 0.5cos(7rt/TU)),  0 < t < Tw, 

w(t)  = < 1,  Tw  < t < T — Tw,  (5.19) 

[ 0.5  — 0.5cos[7r(t  — T)/Tw],  T — Tw  < t < T, 

with  Tw  = T/ 10.  In  the  following  simulations,  we  use  N = 128,  7 = n x 105,  the 

signal  bandwidth  Bs  = 'yT/n,  and  the  sampling  frequency  fs  = 8 Bs.  T is  chosen  in 

such  a way  that  T = Ntt / 87  = 12.6  ms,  Bs  = 1.26  kHz,  fc  = 2 Bs  = 2.53  kHz, 

fs  = 10.12  kHz,  Ts  = 98.82  /is.  For  the  received  signal,  r = 40.5TS  and  /3  = 0.5. 
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Figure  5.2:  Application  of  PEARS  to  a simulated  signal,  (a)  Transmitted  signal,  (b) 
Received  signal  in  the  absence  of  noise,  (c)  Cost  functions  in  the  absence  of  noise, 
(d)  Comparison  of  RMSEs  of  the  time  delay  estimates  and  the  CRB. 


The  SNR  is  defined  as  101og10  where  Es  is  the  energy  of  the  signal  s(nTs).  The 
RMSEs  are  obtained  through  500  independent  Monte  Carlo  trials. 

The  waveforms  of  the  transmitted  signal  and  the  noise  free  received  signal 
are  compared  in  Figures  5.2(a)  and  5.2(b).  The  normalized  cost  functions  for  the 
nonnegative  /3  assumption  (solid  line,  corresponding  to  (5.13))  and  the  complex- 
valued P assumption  (dashed  line,  corresponding  to  (5.17))  are  compared  with  each 
other  in  Figure  5.2(c).  Note  that  the  cost  function  obtained  by  assuming  the  true 
nonnegative  /?  to  be  complex-valued  is  approximately  the  envelope  of  the  true  cost 
function  C\(t).  Note  also  that  we  can  obtain  a more  accurate  estimate  by  maximizing 
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c\(t)  than  maximizing  C2(r)  because  the  former  has  much  sharper  peaks  than  the 
latter.  (In  the  absence  of  noise,  these  two  cost  functions  share  the  same  global 
maximum.)  In  Figure  5.2(d),  the  RMSEs  of  the  estimates  are  compared  with  the 
CRB,  which  is  the  best  performance  bound  of  an  unbiased  estimator.  (The  CRBs 
are  derived  in  Appendix  C.)  Two  estimation  algorithms  are  used  here.  One  is  our  fast 
two-stage  estimation  algorithm  PEARS  and  the  other  is  the  direct  brute  force  search 
method,  referred  to  as  “Search”  in  Figure  5.2(d).  The  Search  method  is  performed  in 
two  stages  as  well.  First,  we  obtain  the  initial  estimate  tq  by  maximizing  C2(r),  and 
next,  we  refine  it  based  on  c^r).  Both  search  procedures  are  performed  by  using  the 
fmin  function  in  Matlab.  It  can  be  seen  from  Figure  5.2(d)  that  both  algorithms 
can  approach  the  CRB  as  the  SNR  increases.  Note  also  that  the  performance  of  the 
Search  method  is  slightly  better  than  PEARS  when  the  SNR  is  low  to  moderate. 
However,  PEARS  is  computationally  much  more  efficient  than  Search. 

5.5  Experimental  Acoustic  Proximity  Ranging  System 

In  this  section,  we  first  describe  the  experimental  setup  used  to  demonstrate 
the  PEARS  algorithm. 

5.5.1  Experimental  Setup 

Experiments  have  been  conducted  by  using  the  commercially  available  Pana- 
sonic ceramic  transducers  (EFRRSB40K5  and  EFRTSB40K5)  with  the  3 dB  beamwidth 
of  about  30°.  Table  1.1  lists  the  parameters  of  the  transducers  used  in  our  exper- 
iments [86].  The  block  diagram  of  the  experimental  setup  for  the  acoustic  ranging 
system  is  shown  in  Figure  5.3.  The  centers  of  the  two  transducers  are  17  mm  apart 
and  they  are  mounted  above  an  aluminum  plate  with  a very  smooth  surface,  which 
is  used  to  simulate  the  planar  acoustically  hard  boundary.  In  order  to  reduce  the 
secondary  reflection,  the  sound  absorption  material  (plastic  foam)  is  used  to  cover 
the  region  around  the  transducers. 
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Figure  5.3:  Experimental  setup  for  the  acoustic  ranging  system. 

A laser  displacement  sensor  (Micro-epsilon  ILP-2000-10)  is  placed  close  to  the 
transducers  with  its  beam  pointing  vertically  to  the  surface  of  the  aluminum  plate. 
This  laser  displacement  sensor  works  in  a range  of  58  ± 5 mm  with  an  accuracy  of 
0.5  /rm  based  on  a triangulation  principle  [78].  Both  the  acoustic  transducers  and 
the  laser  displacement  sensors  are  fixed  on  a 3-D  traverse  controlled  by  a personal 
computer  (PC),  which  is  used  to  adjust  the  position  of  the  transducers  related  to  the 
boundary  with  an  accuracy  of  0.4  ^m.  The  aluminum  plate  is  fixed  horizontally  on  a 
shaker  controlled  by  a function  generator,  which  can  apply  controlled  vertical  vibra- 
tions to  the  object  supported.  The  excitation  signal  for  the  transmitter  is  generated 
by  an  arbitrary  waveform  generator.  A 16-bit  ADC  is  used  to  sample  (fa  = 196, 608 
Hz)  the  output  signal  from  the  laser  displacement  sensor,  the  received  signal  from 
the  receiver,  and  the  excitation  signal  for  the  transmitter.  The  recorded  data  are 
transferred  to  a host  PC  and  the  PEARS  algorithm  is  applied  to  the  recorded  data 
to  determine  the  distance  estimates. 
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The  sound  speed  is  obtained  by  [2] 

c = y/yRTr,  (5.20) 


where  7 = 1.4,  R = 287  J/(kg  • K)  for  the  ideal  gas,  and  Tr  = (273.16  + 21)  K in  our 
experiment.  Thus  we  have  c = 344  m/s. 

5.5.2  Excitation  Waveform  Design 


The  PEARS  algorithm  can  be  applied  to  arbitrary  transmitted  waveforms 
as  long  as  the  waveform  is  periodic.  In  this  experiment,  to  increase  the  parameter 
update  rate  and  decrease  the  signal  peak-to-average  ratio,  the  special  designed  CW 
waveform  proposed  in  Chapter  4 is  adopted.  A sum  of  eight  sinusoids  with  different 
frequencies  but  the  same  amplitude  is  used  to  continuously  excite  the  transmitter. 
The  frequencies  of  the  sinusoids  are  evenly  spaced  in  the  frequency  domain  as  = 
fo  + kAf,  k = 0, 1, ..,  7,  where  /0  and  A / are  chosen  to  satisfy 


A/< 


1 

5 

'max 


(5.21) 


where  rmax  denotes  the  maximum  possible  time  delay  in  a practical  measurement 
scenario,  and 


fo 

A/ 


= Q, 


(5.22) 


with  q being  an  integer,  and 


(5.23) 


with  N being  an  integer  as  well.  Thus,  the  excitation  signal  is  a periodic  signal 
with  the  period  T — 1/A / = rmax,  and  in  each  period,  the  number  of  samples  is 
guaranteed  to  be  the  integer  N. 

According  to  the  requirements  from  (5.21)  to  (5.23),  considering  the  fact  that 
the  maximum  distance  is  about  100  mm  for  our  system,  as  well  as  taking  into  account 
the  resonant  frequency  range  of  the  Panasonic  ceramic  transducers,  we  set  fo  = 
33,792  Hz  and  A / = 1,536  Hz,  respectively.  Thus,  the  period  of  CW  excitation 
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signal  is  0.651  ms  and  the  number  of  samples  for  each  period  is  N — 128.  Further,  to 
reduce  the  peak-to-average  ratio  of  the  transmitted  signal,  the  initial  phase  for  each 
sinusoid  is  chosen  as  {0, 0, 0,  n,  0, 7r,  7r},  respectively. 

5.5.3  Experimental  Results 

We  first  illustrate  the  measured  waveforms.  Then  we  describe  the  experimen- 
tal measurement  results  by  using  PEARS. 

Figure  5.4  shows  the  recorded  signals  where  four  periods  (i.e.,  512  samples) 
are  illustrated.  Figure  5.4(a)  is  the  excitation  signal  for  the  transmitter  based  on 
the  design  in  Section  5.5.2.  Figure  5.4(b)  is  the  transmitted  waveform  when  the 
transmitter  and  the  receiver  are  face  to  face.  The  distance  between  the  transducers 
has  been  measured  precisely.  Thus  this  waveform  is  taken  as  the  known  transmit- 
ted waveform  corresponding  to  the  signal  in  Figure  5.1(b),  and  one  period  of  the 
signal  in  Figure  5.4(b)  (i.e,  128  samples)  is  viewed  as  the  transmitted  signal  used 
in  data  model  in  (5.2).  Figure  5.4(c)  shows  the  crosstalk  between  the  receiver  and 
the  transmitter,  which  is  obtained  when  the  two  transducers  are  closely  spaced  with 
their  beams  pointing  in  the  same  direction  in  the  absence  of  a reflector.  When  there 
is  a reflector,  the  received  signal  is  a superposition  of  the  true  reflected  signal  from 
the  reflector  and  the  crosstalk.  Since  the  crosstalk  is  unchanged  and  independent  of 
the  reflected  signal,  it  can  be  recorded  and  stored  in  the  system  in  advance  and  sub- 
tracted out  from  each  received  signal  thereafter.  Figure  5.4(d)  is  the  received  signal 
when  the  transducers  face  the  boundary.  Figure  5.4(e)  is  the  signal  in  Figure  5.4(d) 
after  crosstalk  subtraction,  and  Figure  5.4(f)  shows  the  discrete  Fourier  transform 
(magnitude)  of  the  signal  in  Figure  5.4(e). 

Consider  now  the  two  cases  of  stationary  level  measurements  for  a fixed  bound- 
ary. First,  we  change  the  distance  between  the  transducers  to  the  boundary  at  46 
different  levels  with  the  distance  difference  between  two  consecutive  levels  being  2 
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mm.  Thus,  the  entire  measurement  distance  between  the  surface  of  the  transduc- 
ers and  the  boundary  covers  from  10  mm  to  100  mm.  The  observation  time  for 
the  acoustic  ranging  system  at  each  level  is  0.5  second,  during  which  we  can  obtain 
1/(2T)  — 768  distance  estimates.  (T  is  the  period  of  the  acoustic  signal  and  is  equal 
to  0.651  ms  in  our  experiment.)  The  averages  of  the  estimates  at  each  level  are  cal- 
culated. Figure  5.5(a)  shows  the  measurement  results.  The  RMSE  of  the  estimates 
obtained  by  the  acoustic  ranging  system  is  0.25  mm  with  respect  to  the  distances 
determined  by  the  traverse.  Note  that  this  error  is  very  small  compared  with  the 
wavelength  of  8.6  mm  of  the  transmitted  signal  for  our  system.  Note  also  that  the 
outputs  of  the  laser  displacement  sensor  are  saturated  at  both  the  lower  and  higher 
distances  due  to  its  limited  working  range  of  10  mm. 

Second,  to  further  compare  the  measurement  accuracy  between  the  acoustic 
ranging  system  and  the  laser  displacement  sensor,  we  make  the  fine  measurements 
within  the  working  range  of  the  laser  displacement  sensor.  For  this  case,  the  distance 
difference  between  two  consecutive  levels  is  reduced  to  0.25  mm  and  41  different 
levels  are  measured.  Thus,  the  entire  measurement  distance  range  covers  only  about 
10  mm.  The  averages  of  the  measurements  are  shown  in  Figure  5.5(b).  The  RMSEs 
of  the  acoustic  measurements  and  the  laser  measurements  are  0.15  mm  and  0.0036 
mm,  respectively. 

Consider  next  the  two  cases  of  continuous  measurements  for  a fixed  bound- 
ary. In  these  cases,  the  position  of  the  transducers  vertically  increases  or  decreases 
with  a constant  speed  within  a certain  range,  and  the  signals  from  the  acoustic  rang- 
ing system  and  the  laser  displacement  sensor  are  recorded  simultaneously  during 
such  movements.  The  observation  time  is  one  second,  during  which  we  can  obtain 
1/^—1536  estimates.  Figure  5.6(a)  shows  the  measurement  results  for  the  first  case 
when  the  distance  decreases  from  110  mm  to  10  mm  in  one  second.  The  true  distance 
(green  solid  line)  in  this  figure  is  obtained  by  first-order  curve  fitting  according  to 
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the  outputs  of  the  laser  displacement  sensor  in  its  working  range.  The  RMSE  of  the 
acoustic  measurements  is  0.27  mm.  Figure  5.6(b)  shows  the  measurement  results  for 
the  second  case  when  the  measurement  is  performed  mainly  in  the  working  range  of 
the  laser  displacement  sensor.  The  RMSE  of  the  acoustic  measurements  with  respect 
to  the  laser  measurements  is  0.12  mm.  Note  that  in  the  above  processing,  only  the 
first  estimate  for  each  case  is  obtained  by  using  the  two-stage  PEARS  algorithm. 
The  estimates  that  follow  are  obtained  by  using  the  simplified  method  proposed  in 
Section  5.3  (i.e,  the  current  refined  estimate  is  used  as  the  initial  estimate  for  the 
next  estimation  with  the  second  stage  of  PEARS). 

It  can  be  noted  from  Figures  5.5(b)  and  5.6(b)  that  the  measurements  obtained 
by  the  acoustic  ranging  system  slightly  fluctuate  around  the  true  distances.  This 
error,  in  fact,  is  caused  mainly  by  the  secondary  reflection  of  the  acoustic  wave  via  the 
surface  of  the  transducers,  which  results  in  a small  portion  of  a standing  wave  pattern 
being  formed  between  the  surface  of  the  transducers  and  the  boundary.  However, 
this  error  can  be  effectively  removed  through  calibration,  for  example,  by  using  a 
calibration  curve  obtained  via  a polynomial  fitting  based  on  the  real  measurement 
bias. 

Finally,  PEARS  is  applied  to  process  the  experimental  data  recorded  when 
the  boundary  is  nonstationary  (i.e.,  when  the  shaker  is  vibrating).  The  observation 
time  for  each  case  is  one  second.  The  measurements  obtained  by  the  acoustic  rang- 
ing system  have  been  calibrated  based  on  the  calibration  curve  obtained  from  the 
measurement  results  in  Figure  5.6(b). 

In  our  experiments,  three  different  driving  signals  are  applied  to  the  shaker. 
For  the  first  case,  the  amplitude  and  the  frequency  of  the  driving  signal  are  1 V and 
2 Hz,  respectively.  Figure  5.7(a)  shows  the  measurement  results.  The  RMSE  of  the 
acoustic  measurements  with  respect  to  the  laser  measurements  is  0.019  mm.  For 
the  second  case,  the  amplitude  and  the  frequency  of  the  driving  signal  to  the  shaker 
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are  1.5  V and  10  Hz,  respectively.  Figure  5.7(b)  shows  the  measurement  results, 
and  the  corresponding  RMSE  is  0.031  mm.  For  the  third  case,  the  amplitude  and 
the  frequency  of  the  driving  signal  to  the  shaker  are  0.5  V and  20  Hz,  respectively. 
Figure  5.7(c)  shows  the  measurement  results,  and  the  corresponding  RMSE  is  0.0096 
mm.  As  can  be  seen  from  Figure  5.7,  the  results  obtained  by  the  acoustic  system 
and  the  laser  system  are  well  matched  and  no  obvious  discrepancies  between  the  two 
measurement  methods  are  observed. 

5.6  Summary 

In  this  chapter,  a fast  two-stage  time  delay  estimation  algorithm,  which  can 
be  applied  to  arbitrary  waveforms  transmitted  and  is  referred  to  as  PEARS,  is  pro- 
posed. Numerical  results  demonstrate  that  the  RMSEs  of  the  estimates  obtained  via 
PEARS  can  achieve  the  CRB  as  the  SNR  increases.  Experimental  results  based  on 
the  commercially  available  transducers  show  that  our  measurement  method  can  be 
used  to  obtain  a high  measurement  accuracy  for  a ranging  distance  of  less  than  100 
mm  with  the  parameter  update  rate  up  to  around  1.5  kHz. 
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Figure  5.4:  Recorded  waveforms  for  the  experiment,  (a)  The  waveform  of  the  excita- 
tion signal  applied  to  the  transmitter,  (b)  The  transmitted  waveform,  (c)  Crosstalk, 
(d)  The  received  signal  when  the  transducers  face  the  boundary,  (e)  The  signal  in 
(d)  after  the  crosstalk  subtraction,  (f)  Discrete  Fourier  transform  of  the  signal  in  (e). 
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Figure  5.5:  Stationary  level  measurements  for  a fixed  boundary,  (a)  Ranging  distance 
covers  from  10  mm  to  100  mm.  The  distance  difference  between  two  consecutive  levels 
is  2 mm.  (b)  Fine  measurements  within  the  working  range  of  the  laser  displacement 
sensor.  The  distance  difference  between  two  consecutive  levels  is  0.25  mm. 
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(b) 

Figure  5.6:  Continuous  measurements  for  a fixed  boundary,  (a)  Measurement  results 
when  the  distance  decreases  from  110  mm  to  10  mm  with  a constant  speed,  (b) 
Measurement  results  when  the  distance  increases  from  45  mm  to  55  mm  with  a 
constant  speed. 
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Figure  5.7:  Comparison  of  the  displacement  measurements  obtained  by  the  acoustic 
ranging  system  and  the  laser  displacement  sensor.  Measurement  results  when  the 
amplitude  and  frequency  of  the  input  signal  to  the  shaker  are  (a)  1 V and  2 Hz,  (b) 
1.5  V and  10  Hz,  and  (c)  0.5  V and  20  Hz,  respectively.  (Note  that  the  dash-dot  line 
and  solid  line  almost  coincide  with  each  other.) 


CHAPTER  6 

MULTI-PEARS  ALGORITHM 

6.1  Introduction 

In  the  experiment  section  of  Chapter  5,  we  notice  that  when  ultrasonic  sensors 
face  a sound-hard  reflective  surface  with  a very  small  distance,  the  reflected  sound 
wave  can  bounce  back  and  forth  several  times  between  the  sensors  and  the  reflection 
surface  before  decaying  to  zero,  which  results  in  unwanted  strong  and  overlapping 
secondary  echoes  in  the  received  signal.  The  time  delays  for  the  secondary  echoes  are 
approximately  integer  multiples  of  the  time  delay  of  the  first  echo.  The  existence  of 
secondary  echoes  can  degrade  the  performances  of  our  former  developed  algorithms 
because  of  the  model  mismatch.  Although  some  simple  schemes  can  be  adopted  to 
deal  with  the  secondary  echoes,  such  as  reducing  the  secondary  reflection  area  by 
covering  it  with  the  sound  absorption  materials  and  using  the  calibration  to  remove 
the  estimation  bias  caused  by  secondary  echoes,  the  new  measurement  techniques 
which  take  the  secondary  echoes  into  account  need  to  be  investigated. 

The  traditional  matched  filter  based  methods  cannot  resolve  two  echoes  with 
a time  spacing  less  than  the  reciprocal  of  the  signal  bandwidth  [109].  Hence,  for  most 
of  the  very  short  distance  measurement  scenarios,  the  matched  filter  based  algorithms 
tend  to  fail  or  suffer  from  severe  performance  degradations  due  to  their  poor  reso- 
lutions [34].  Most  of  the  existing  super- resolution  time  delay  estimation  approaches 
(see  [110,  111]  and  the  references  therein)  are  developed  for  general  purposes.  They 
do  not  exploit  the  a priori  information  of  the  integer  multiple  time  delays  and  the 
nonnegative  amplitude  due  to  the  acoustically  hard  reflections. 
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In  this  chapter,  we  focus  on  acoustic  proximity  ranging  in  the  presence  of 
secondary  echoes.  A new  time  delay  estimation  method  is  proposed  for  the  joint 
proximity  ranging  and  secondary  echo  mitigation.  We  establish  a flexible  data  model 
that  takes  into  account  the  acoustically  hard  reflection  as  well  as  small  random  time 
delay  variations  of  the  secondary  echoes  around  the  integer  multiples  of  the  time  delay 
of  the  first  echo.  The  new  time  delay  estimator  is  based  on  the  NLS  fitting  criterion. 
However,  it  is  difficult  to  directly  optimize  the  highly  nonlinear  NLS  cost  function.  We 
present  a novel  computationally  efficient  time  delay  estimation  algorithm,  referred  to 
as  Multi-PEARS  (Multi-echo  Parameter  Estimation  for  Acoustic  Ranging  Systems), 
to  optimize  the  NLS  cost  function. 

6.2  Data  Model  and  Problem  Formulation 

As  shown  in  Figure  1.1(b),  a two-transducer  ranging  system  is  used  to  measure 
the  distance  between  the  transducers  and  an  acoustically  hard  boundary.  We  assume 
that  the  propagation  medium  is  homogeneous  and  nondispersive  and  the  common 
region  covered  by  both  the  transmitter  and  the  receiver  is  a smooth  planar  surface. 
The  transmitted  signal  is  emitted  from  the  transmitter  toward  the  reflection  bound- 
ary. The  directly  reflected  signal  received  by  the  receiver  is  called  the  first  echo. 
Since  the  distance  between  the  sensors  and  the  boundary  is  very  short,  the  secondary 
reflections  between  the  sensors  and  the  ranging  surface  cannot  be  ignored.  The  time 
delays  of  the  secondary  echoes  are  approximately  integer  times  of  the  time  delay  of 
the  first  echo.  Hence  the  received  signal  can  be  modeled  as 

L 

y(t)  = ^2ais(t- It  - At) +e(t),  0 < t < T,  (6.1) 

i=i 

where  s(f),  0 < t < T,  represents  the  known  (or  measured)  real- valued  transmit- 
ted signal  and  T is  the  signal  duration;  y(t ) denotes  the  real-valued  received  signal; 
{«/}(= i denote  the  amplitudes  of  the  received  signal  and  is  assumed  to  be  nonneg- 
ative due  to  the  acoustically  hard  reflection;  r denotes  the  time  delay  corresponding 


78 


to  the  first  echo,  and  the  time  delays  for  the  secondary  echoes  are  It  + A;,  l = 2, L, 
with  Ai  = 0 and  {A i}[=2  being  small  perturbations  around  the  integer  multiples  of 
r;  e(t)  is  the  real-valued  noise,  which  is  modeled  as  a zero-mean  Gaussian  random 
process. 

The  received  signal  after  sampling  and  A/D  conversion  has  the  form: 

L 

y{nTs)  = ^2ais(nTs  - It  - At)  + e(nTs),  n = 0,1,  ...,N  - 1,  (6.2) 

j=i 

where  Ts  denotes  the  sampling  period,  which  is  equal  to  the  reciprocal  of  the  sampling 
frequency  fs. 

Our  problem  of  interest  herein  is  to  estimate  r from  {y(nTs)}^T01  when  s(t), 
0 < t < T,  or  more  practically  {s(nTs)}„~Q  is  known. 

Instead  of  working  on  this  problem  in  the  time  domain,  the  frequency  domain 
is  preferred.  This  is  due  to  the  following  two  reasons.  First,  for  the  time  domain 
methods,  we  could  either  restrict  r to  be  multiples  of  Ts  or  resort  to  interpolation  if 
only  {s(nTs)}^701  is  known  [111].  This  inconvenience  can  be  avoided  by  transforming 
the  problem  into  the  frequency  domain,  where  r can  take  on  a continuum  of  values. 
Second,  for  the  time  domain  data  model  in  (6.2),  { on , A i}fL2  for  each  secondary  echo 
must  be  dealt  with  separately.  However,  by  transforming  (6.2)  into  the  frequency 
domain  and  using  the  facts  that  {A /}/L2  are  small  values,  a more  concise  data  model 
can  be  obtained. 

Let  Y(m),  S(m),  and  E(m),  where  m = -A/2,  -A/2  + 1, ...,  A/2  - 1,  denote 
the  DFTs  of  y(nTs),  s(nTs),  and  e(nTs),  respectively.  Provided  that  the  aliasing  due 
to  sampling  is  negligible,  Y (m)  can  be  written  as 

L 

y(m)  = a/S(m)e-^(/T+A')m  + E{m).  (6.3) 

(=i 

Since  both  the  transmitted  signal  s(t)  and  the  received  signal  y(t)  are  real-valued, 
their  Fourier  transforms  are  conjugate  symmetric.  Most  of  the  currently  available 
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ultrasonic  transducers  are  resonant  devices  (i.e.,  underdamped  second-order  dynamic 
systems)  [17].  Hence  the  transmitted  signal  s(t)  is  narrowband  in  general.  In  the 
negative  frequency  domain,  i.e.,  for  m — —N/2,  —N/2  + 1, ...,  0,  since  s(t ) is  narrow- 
band,  S(m)  mainly  occupies  a few  frequency  bins  in  the  range  [m0  — Am,  m0  + Am], 
where  m0  is  the  peak  location  of  ^(m)!  and  2 Am  « BSNTS  with  Bs  being  the  signal 
bandwidth.  Based  on  the  above  assumption  and  for  very  small  {A i}£_2,  we  have  the 
following  approximations  for  the  signal  term  in  (6.3): 

aiS(m)e-j^ 7(,T+A')m  = a/5(m)e'J^A,(m"mo)e-^A,m°e~J^iTm 

« a/5(m)e-^A,moe-*,Tm 
= ^S{m)e-j^lTm,  (6.4) 


where  fa  = aie_J7^A'm°,  / = 1, 2, ...,  L.  Note  that  fa  = au  is  nonnegative  and  {$}/L2 
are  complex- valued.  In  the  second  step  of  deriving  (6.4),  we  have  assumed  that 
S(m)  « 0 for  m [mo  — Am,  mo  + Am]  and  {A are  so  small  that  « 0 

(i.e.,  'kAiBs  « 0),  and  hence  « 1 for  m G [m0— Am,  mo+Am].  There- 

fore the  small  time  delay  variations  of  the  secondary  echoes  have  been  transformed 
to  the  phases  of  their  complex  amplitudes. 

Based  on  (6.3)  and  (6.4)  and  by  exploiting  the  conjugate  symmetry  property 
of  Y(m),  S(m)  and  E(m),  the  unknown  parameters  can  be  obtained  by  minimizing 
the  following  NLS  criterion: 


A({/3,r})  = 


E p2(™) 


L 

Y(m)  - S(m)^2fa 


2 


. 2kIt 
J NT a 


(6.5) 


where  fa  = [ fa  fa 
and  P{-N/ 2)  = P( 0)  = 
Define 


m=-N/2 


1=1 


fa  ]T  with  (-)r  denoting  the  transpose,  {P(m)  = l}mL_N/2+i’ 


i 

Vi’ 


x = [ P(-N/2)Y(-N/2)  P(-N/2  + l)Y(-N/2  + l)  •••  P(0)T(0)  ]T,  (6.6) 
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where  K = N/2  + 1 is  the  dimension  of  the  vector  x.  Let  W be  a K x K diagonal 
matrix: 

W = diag{P(— iV/2)S(— iV/2),  P(—N/2  + l)S(-N/2  + 1),  • • • , P(0)S(0)}.  (6.7) 

Let 

b (It)  = [e-J^(-N/2)  e-jffi(-N/2+1)  ■■■  1]T 

_ [ e~j2nf0lr  e~j2irfilT  , . . e~j2nfK-ilT  jT^ 

where 

fk  = ~2Ts  + NTS'  = 

Define 

B(r)  = [b(r)  b(2r)  •••  b(Lr)  ], 

and 

D(r)  = WB(r). 

Then  (6.5)  can  be  expressed  as 

£i({/3,r})  = ||x-  D(r)^||2, 

where  ||  • ||  denotes  the  Euclidean  norm.  For  the  white  noise  case,  the  above  NLS 
approach  is  the  same  as  the  ML  method.  When  the  noise  is  colored,  the  NLS  approach 
can  still  give  very  accurate  parameter  estimates  [65]. 

We  remark  that  the  time  delay  estimation  problem  considered  herein  is  similar 
to  the  well-known  harmonic  sinusoidal  parameter  estimation  problem  [62,  88,  114]. 
However,  the  sinusoidal  signals  in  our  problem  are  weighted  by  the  known  signal 
spectrum  and  the  amplitude  for  the  first  echo  is  nonnegative.  Since  most  of  the 
existing  harmonic  retrieval  methods  are  designed  for  signals  with  complex-valued 
amplitudes  and  flat  band-limited  spectra,  they  are  not  directly  applicable  to  our 
problem  of  interest. 


(6.8) 

(6.9) 

(6.10) 

(6.11) 

(6.12) 
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Note  that  minimizing  (6.12)  is  a highly  nonlinear  optimization  problem.  We 
present  below  an  efficient  optimization  method,  referred  to  as  the  Multi-PEARS 
algorithm,  to  minimize  (6.12)  by  taking  into  account  both  nonnegative  P\  and  the 
weighted  harmonic  structure  of  the  data  model. 

6.3  The  Multi-PEARS  Algorithm 

The  Multi-PEARS  algorithm  optimizes  the  NLS  cost  function  iteratively. 
First,  we  assume  all  of  the  amplitudes  {Pi}^  are  complex-valued  and  obtain  an 
initial  estimate  by  minimizing  (6.12)  via  an  FFT-based  method.  Next,  based  on  the 
initial  estimates,  we  subtract  out  the  components  of  the  secondary  echoes  and  keep 
only  the  component  of  the  first  echo.  Then,  the  single-echo  based  PEARS  algorithm 
proposed  in  Chapter  5 is  used  to  refine  the  initial  time  delay  estimate.  The  refining 
step  is  iterated  several  times. 

6.3.1  Initialization  Stage 

Assume  that  {Pi}i=x  are  complex- valued,  it  can  be  shown  that  the  unknown 
parameters  minimizing  (6.12)  can  be  determined  by  [99]: 

f0  = argmaxco(r),  (6.13) 

T 

where 

c0(r)  = xhD(t  ) [DH(r  )D(r  )]-1DH(r  )x,  (6.14) 

and 

00  = |D«(T)D(r)]-1D'I(r)x|„^,  (6.15) 

with  (• )H  denoting  the  conjugate  transpose.  According  to  the  parsimony  principle 
[97],  if  P\  is  indeed  nonnegative,  the  estimate  obtained  by  maximizing  c0(r)  should 
be  less  accurate  than  that  obtained  by  minimizing  the  true  cost  function  of  (6.12). 
Note  that  (6.13)  is  a one-dimensional  search  problem  and  co(r)  usually  has  many 
local  maxima.  Hence,  maximizing  c0(r)  requires  the  search  over  a very  fine  grid.  By 
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taking  advantage  of  the  harmonic  structure  of  the  data  model,  we  use  zero-padding 
FFT  to  perform  such  a search  computationally  efficiently. 

Note  that  D//(r)D(r)  in  (14)  can  be  expressed  as 


D*(t)D(t)  = 


do  d\  ...  di- i 

dt  ; 

: •••  •••  a, 

d*L_i  ...  d\  d0  _ 


LxL 


(6.16) 


where  (•)*  denotes  the  complex  conjugate  and 


K- 1 


dt  = Y,\w(k)\2e~j2nfklT 


k= 0 


N/2 


e^^\W(k)\2e-juk 


k= 0 


, l = 0,1, L — 1, 


(6.17) 


2irlr 

U~NTl 


with  W{k)  denoting  the  /cth  diagonal  element  of  the  diagonal  matrix  W in  (6.7),  and 
xHD(r)  can  be  expressed  as 


xhD(t)  = [ ai  a2  • • • aL  ], 


(6.18) 


where 


K- 1 


ai  = Y^[x*(k)W(k)}e~j27rfklT 


k= 0 


N/2 


= eJ'^^[x*(A:)fF(fc)]e-J'a'fc 


fc=0 


(6.19) 


27 r/r 

u~nt; 


with  x(/c)  denoting  the  fcth  element  of  x in  (6.6). 

It  can  be  seen  from  (6.17)  and  (6.19)  that  {di}fZo  and  {a;}^  can  be  ob- 
tained as  follows.  First,  the  DTFT,  which  can  be  efficiently  implemented  by  using 
zero-padding  FFT,  is  applied  to  the  sequences  {|IT(fc)|2}^o  and  {x* (k)W (k)}^l 20, 
respectively.  Next,  the  linear  phase  term  ^ is  multiplied  to  the  outputs  of  DTFT, 
and  finally,  the  sampling  is  performed  at  the  frequency  points  u — ^r. 
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Define  the  following  two  vectors  which  are  obtained  by  multiplying  the  phase 
term  e 4^  to  the  FFTs  of  {|fT(fc)|2}^g  and  {x*  (k)W  (k)}^L q,  respectively,  as 


q = [ q{ 0) 

<7(1)  • 

••  q(N- 

- 1)  f, 

(6.20) 

r = [ r(0) 

r(l) 

■■  r(N  - 

- 1)  ]T, 

(6.21) 

where  N denotes  the  data  length  after  zero-padding.  For  N sufficient  large  and  for  any 
given  r and  l,  the  corresponding  d;  and  a*  can  be  approximately  obtained  by  picking 
the  elements  of  q and  r at  the  index  ( Nh)/(NTS ) , where  |yj  denotes  rounding  to 
the  nearest  integer.  Once  {dJ/Fg1  and  {a*}^  are  in  hand,  Dh(t)D(t)  and  xHD(r) 
can  be  constructed  and  substituted  into  (6.14)  to  obtain  the  cost  function  Co(r)  at 
the  r.  The  calculation  steps  for  the  initial  estimation  can  be  summarized  as  follows: 


Step  (1):  Zero-pad  {\W  {k)\2}^=Q  and  {x*  (k)W  (k)}%=0  to  length  of  N and  perform  Ap- 
point FFTs.  Compensate  out  the  linear  phase  term  for  the  FFTs  and 
obtain  the  data  vectors  q and  r,  respectively. 


Step  (2):  For  a given  r,  obtain  di  and  a/  as 


di  = q{  [(NtI)/(NTs)\ ) , l = 0, 1,  • • • , L — 1, 


(6.22) 


and 


ai  = r([(NTl)/(NT,)  J),  / = 1,2,  ...,L. 


(6.23) 


Step  (3):  Construct  D//(r)D(r)  and  xHD(r)  by  using  {d;}^1  and  {fo}^  obtained  in 
Step  (2).  Compute  cq(t)  for  the  given  r based  on  (6.14). 


Step  (4):  According  to  the  maximum  peak  location  of  the  computed  c0(r)  in  Step  (3), 
obtain  the  initial  estimate  f0  of  r.  Once  f0  is  obtained,  (30  can  be  readily 
calculated  by  using  (6.15). 
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In  practical  applications,  we  have  used  N = 4IV,  which  is  sufficiently  accurate  for  the 
initial  estimation  in  our  examples.  Note  also  that  in  the  above  search  method,  we  need 
to  compute  the  FFTs  of  {\W[k)\2}^Q  and  {x* (k)W 20  only  once.  Thereafter, 
the  look-up  table  method  is  used  to  find  {d/}^1  and  for  a given  r.  Thus, 

this  method  is  computationally  more  efficient  than  the  brute  force  search  method, 
and  can  be  easily  implemented  by  currently  available  dedicated  FFT  chips. 

6.3.2  Refining  Stage 

Assume  {$i}[=2  and  f are  given,  we  subtract  out  the  secondary  echoes  to 

obtain 

L 

xi  = x — ^AWb(lf).  (6.24) 

1=2 

Then  (6.12)  becomes 

A({A,t})  = ||Xl  - AWb(r)||2,  (6.25) 

where  Pi  is  assumed  to  be  nonnegative  due  to  the  acoustically  hard  reflection.  We 
have  developed  an  efficient  optimization  algorithm,  referred  to  as  PEARS,  to  solve 
the  time  delay  estimation  problem  for  this  single  echo  case  [67].  Note  that  PEARS 
is  a two-stage  algorithm  which  first  obtains  an  initial  estimate  based  on  a smooth 
cost  function,  and  then  refines  the  initial  estimate  based  on  the  true  cost  function. 
Herein,  since  we  have  already  obtained  the  initial  estimate  of  r,  we  only  need  to  use 
the  second  stage  of  PEARS  to  refine  it.  Note  that  iteration  can  improve  the  accuracy 
of  f . The  refined  estimate  f can  be  used  to  redetermine  (3  via  (6.15).  Then  the  refined 
f and  (3  are  used  to  redetermine  xx  by  using  (6.24).  Based  on  the  updated  xi,  a more 
accurate  estimate  f of  r can  be  obtained  by  using  PEARS.  With  the  above  simple 
preparations,  we  now  present  the  steps  of  the  Multi-PEARS  algorithm. 

Step  1:  Obtain  the  initial  estimates  f0  and  (30  by  using  the  FFT-based  method  proposed 
in  Section  6.3.1.  Let  t — t0  and  (3  = {30. 


Step  2:  Obtain  Xi  by  using  (6.24). 
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Step  3:  Obtain  the  refined  f by  using  the  second  stage  of  PEARS  as  follows: 


Substep  3.1:  Obtain  the  measurement  for  each  frequency  point  as  Xk  — W*(k)x\{k), 
k = 0,1, ...,  K — 1,  where  xx(k ) denotes  the  /cth  element  of  xx. 

Substep  3.2:  Obtain  the  phase  for  each  measurement  by  h = ~ arg(xfc),  k = 0, 1, K- 
1,  where  arg(:r)  denotes  the  phase  of  x. 


Substep  3.3:  Obtain  the  integers  Zk  = 


'Mkl-h 

2tt 


, k — 0, 1, K — 1,  where  fk  is  defined 


in  (6.9)  and  [-J  denotes  rounding  to  the  nearest  integer. 
Substep  3.4:  Obtain  the  refined  t via  f = ^k=Ji x_k} ik ^ k Zk ^ . 

^k=o  2ir\£k\fk 


Step  4:  Determine  refined  (3  via  (6.15)  by  using  the  refined  f. 


Step  5:  Iterate  Steps  2 through  4 and  calculate  the  residue  after  each  iteration  as 
||x  — D(f)/3||2.  If  the  relative  change  of  the  residues  between  two  consecutive 
iterations  is  less  than  a small  value  e,  the  algorithm  stops.  (In  our  examples, 
we  have  used  e = 10-4  to  test  the  convergence  of  the  algorithm.) 

Before  we  proceed,  let  us  remark  on  the  following  two  facts.  First,  Multi- 
PEARS  converges  very  fast  due  to  the  good  initial  estimates.  In  practice,  we  note 
that  the  algorithm  usually  converges  within  3 or  4 iterations.  Second,  the  number  of 
the  echoes  L is  determined  by  the  practical  measurement  environment.  Whenever  L 
is  unknown,  it  can  be  estimated  from  x by  using,  for  instance,  the  generalized  Akaike 
information  criterion  [65]. 

6.4  Numerical  and  Experiment  Results 

In  this  section,  we  first  give  a brief  introduction  to  the  experimental  setup  used 
to  collect  the  measurement  data.  When  then  present  numerical  examples  to  show  the 
excellent  resolution  and  estimation  performance  of  Multi-PEARS,  where  the  RMSEs 
of  the  estimates  obtained  via  Multi-PEARS  are  compared  to  the  corresponding  CRBs. 
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Finally,  the  new  algorithm  is  applied  to  process  the  real  experimental  data  to  show 
its  capability  in  practice. 

6.4.1  Experimental  Setup 

The  experimental  setup  is  similar  to  that  used  in  Section  5.5  except  that 
the  laser  displacement  sensor  is  not  employed  in  this  case.  The  photograph  for  the 
acoustic  ranging  system  is  shown  in  Figure  6.1.  The  geometric  centers  of  the  two 
transducers  are  17  mm  apart  and  they  are  mounted  above  an  aluminum  plate  with 
a smooth  surface,  which  is  used  to  simulate  the  planar  acoustically  hard  boundary. 
The  acoustic  transducers  are  fixed  on  a 3-D  traverse  controlled  by  a PC,  which  is  used 
to  adjust  the  position  of  the  transducers  related  to  the  boundary  with  an  accuracy 
of  0.4  //m. 

The  excitation  signal  for  the  transmitter  is  a pulse  burst  generated  by  an 
arbitrary  waveform  generator.  The  pulse  repetition  frequency  is  192  Hz  and  the 
pulse  width  is  0.15  ms.  Each  pulse  has  been  modulated  by  a carrier  frequency  of 
40  kHz,  which  is  the  resonant  frequency  of  the  transducer.  (In  this  experiment,  the 
periodic  CW  waveform  adopted  in  Chapter  5 is  not  used.  This  is  simply  because  that 
the  change  of  the  secondary  echoes  is  more  apparent  by  using  the  pulse  waveform 
rather  than  the  CW  waveform.  In  the  real  application,  arbitrary  waveforms  can  be 
employed  as  long  as  they  are  periodic.)  A 16-bit  ADC  is  used  to  sample  (/s  = 196, 608 
Hz)  the  excitation  signal  for  the  transmitter  and  the  received  signal  from  the  receiver. 
The  recorded  data  are  transferred  to  a host  PC  and  the  Multi-PEARS  algorithm  is 
applied  to  the  recorded  data  to  determine  the  distance  estimates.  The  sound  speed 
is  obtained  by  (5.20)  and  c = 344  m/s  in  this  experiment. 

Since  the  transmitter  and  the  receiver  are  closely  placed,  there  exists  a small 
amount  of  crosstalk.  Note  that  the  crosstalk  is  unchanged  and  independent  of  the 
reflected  signal,  and  hence  it  can  be  recorded  and  stored  in  the  system  in  advance  and 
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Figure  6.1:  Photograph  of  the  setup  for  the  sensors  and  the  boundary. 


subtracted  out  from  each  received  signal  thereafter.  All  the  received  signals  below 
are  illustrated  after  the  crosstalk  subtraction. 

6.4.2  Numerical  Examples 

To  demonstrate  the  estimation  accuracy  of  the  proposed  Multi-PEARS  algo- 
rithm, we  present  a numerical  example  below  based  on  the  real  transmitted  signal 
as  shown  in  Figure  6.2(a).  This  waveform  is  obtained  when  the  transmitter  and  the 
receiver  are  face  to  face.  The  distance  between  the  transducers  has  been  measured 
precisely.  Thus  this  waveform  is  taken  as  the  known  transmitted  waveform.  Fig- 
ure 6.2(b)  shows  the  discrete  Fourier  transform  (magnitude)  of  the  signal  in  Figure 
6.2(a).  It  can  be  seen  that  the  transmitted  signal  is  narrowband  with  the  bandwidth 
about  4 kHz.  Thus,  the  nominal  minimum  distance  required  to  resolve  the  first  and 
secondary  echoes  by  the  matched  filter  method  is  about  43  mm. 

In  the  simulation,  define  the  reciprocal  of  the  signal  bandwidth  as  Te  = 1/BS  — 
0.25  ms.  For  the  received  signal,  r = 24.5 Ts  = 0.5 Te,  aq  = 0.5,  a-2  = 0.3,  c*3  = 0.1, 
and  A2  = A3  = 0.  Note  that  the  time  delay  between  each  echo  is  only  about  one 
half  of  Te.  Thus,  the  first  echo  and  secondary  echoes  cannot  to  be  resolved  by  using 
the  matched  filter  based  methods  in  this  case.  The  noise  free  received  signal  is  shown 
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Figure  6.2:  Numerical  example  based  on  the  real  transmitted  waveform  in  the  ex- 
periment. (a)  Transmitted  waveform,  (b)  Discrete  Fourier  transform  of  the  signal  in 
(a),  (c)  Received  signal  in  the  absence  of  noise. 

in  Figure  6.2(c).  The  SNR  is  defined  as  101og10  where  Es  is  the  energy  of  the 
signal  s(nTs).  The  RMSEs  are  obtained  through  200  Monte  Carlo  trials. 

In  Figure  6.3,  the  RMSEs  of  the  distance  estimates  [ref  2 with  c = 344  m/s) 
are  compared  with  the  corresponding  CRBs  (the  CRBs  are  derived  in  Appendix 
D).  Two  different  cases  for  the  time  delay  variations  of  the  secondary  echoes  are 
illustrated:  1)  A2  = A3  = 0,  and  2)  A2  = A3  = 0.1Te.  In  case  (1),  the  RMSEs  of  the 
estimates  and  the  corresponding  CRB  are  denoted  as  the  circles  and  the  solid  line, 
respectively.  In  case  (2),  the  RMSEs  of  the  estimates  and  the  corresponding  CRB  are 
denoted  as  the  squares  and  the  dashed  line,  respectively.  It  can  be  seen  that  when 
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Figure  6.3:  Compares  the  RMSEs  of  the  estimates  and  the  corresponding  CRBs 
for  the  distance  (re/ 2 with  c = 344  m/s).  Two  different  cases  for  the  time  delay 
variations  of  the  secondary  echoes  are  illustrated:  1)  A2  = A3  = 0.  The  RMSEs 
and  the  corresponding  CRB  are  denoted  as  the  circles  and  the  solid  line,  respectively. 
2)  A2  = A3  = 0.1Te.  The  RMSEs  and  the  corresponding  CRB  are  denoted  as  the 
squares  and  the  dashed  line,  respectively. 


there  are  no  time  delay  variations  for  the  secondary  echoes  (i.e.,  A2  = A3  = 0),  the 
RMSEs  of  the  estimates  can  achieve  the  corresponding  CRB.  Note  also  the  threshold 
effects  where  the  RMSEs  of  the  estimates  deviate  from  the  CRB  when  the  SNR 
is  less  than  20  dB.  When  there  exist  the  time  delay  variations  for  the  secondary 
echoes,  the  estimates  are  biased  and  the  RMSEs  of  the  estimates  cannot  achieve  the 
corresponding  CRB.  However,  if  the  time  delay  variations  are  very  small  compared 
with  the  reciprocal  of  the  signal  bandwidth,  the  RMSEs  of  the  estimates  can  still 
approach  the  corresponding  CRB,  and  thus  our  method  can  achieve  a high  estimation 
accuracy  in  these  cases  as  well. 

Figures  6.4(a)  to  6.4(c)  illustrate  the  RMSEs  of  the  estimates  compared  with 
the  CRBs  as  the  function  of  time  delay  by  fixing  the  SNR.  The  time  delays  are 
normalized  to  Te  and  we  assume  A2  = A3  = 0.  It  can  be  seen  that  as  the  time 
delay  decreases,  the  corresponding  CRB  increases.  This  is  as  expected  because  the 
smaller  the  time  delay,  the  more  the  echoes  are  overlapped,  and  hence  the  worse  the 
estimation  accuracy.  Noted  also  the  threshold  effects  in  Figures  6.4(a)  and  6.4(b) 
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Figure  6.4:  RMSEs  and  CRBs  vs.  the  normalized  time  delay  at  different  SNRs.  The 
time  delays  are  normalized  to  Te,  which  is  the  reciprocal  of  the  signal  bandwidth,  (a) 
SNR=20  dB,  (b)  SNR=25  dB,  and  (c)  SNR=30  dB. 


where  the  RMSEs  of  estimates  deviate  from  the  CRB  when  the  time  delay  is  less 
than  a certain  value.  In  order  to  still  achieve  the  corresponding  CRB  in  the  smaller 
time  delay  region,  the  higher  SNR  is  required. 

6.4.3  Experimental  Results 

Now  we  apply  Multi-PEARS  to  process  the  experimental  data  collected  by  the 
ranging  system  introduced  previously.  We  increase  the  distance  between  the  trans- 
ducers to  the  boundary  starting  at  29  mm.  The  entire  measurement  range  is  exactly 
10  mm  and  is  covered  in  21  discrete  levels.  The  distance  difference  between  two 
consecutive  levels  is  0.5  mm.  Figures  6.5(a)  to  6.5(d)  show  the  measured  waveforms 
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recorded  at  the  levels  3,  5,  15  and  19,  respectively.  As  we  can  see  that  due  to  the 
presence  of  the  secondary  echoes,  the  shapes  of  the  received  signals  at  different  levels 
are  quite  different. 

Figure  6.6  shows  an  example  where  the  Multi-PEARS  algorithm  is  applied  to 
the  recorded  signal  at  level  19.  We  use  generalized  Akaike  information  criterion  to 
determine  the  number  of  the  echoes  and  obtain  L — 4 for  this  example.  It  can  be 
seen  from  this  figure  that  Multi-PEARS  can  correctly  resolve  the  first  and  secondary 
echoes  even  though  they  are  heavily  overlapped.  Note  also  that  the  reconstructed 
signal  matches  the  original  signal  very  well.  The  relative  error  between  these  two 
signals  is  only  0.7%. 

Figure  6.7(a)  shows  the  measurement  results  for  the  entire  range  with  the 
different  choices  for  L (the  number  of  echoes).  If  the  secondary  echoes  are  not  taken 
into  account  and  L — 1 is  used,  an  accurate  f cannot  be  obtained.  By  using  L — 2, 
the  measurement  results  are  improved.  But  there  are  still  four  levels  where  the 
measurements  are  quite  different  from  the  true  distances.  (The  results  for  L — 3 are 
similar  to  those  of  L — 2.)  By  using  L = 4,  very  good  results  are  obtained  for  all 
distances.  The  zoomed-in  view  of  the  measurement  errors  are  shown  in  Figure  6.7(b). 
Note  that  the  maximum  measurement  error  is  about  0.02  mm  for  L = 4.  The  ratio 
between  this  error  and  the  entire  measurement  range  is  2 x 10-3. 

6.5  Summary 

Due  to  the  presence  of  unwanted  strong  and  closely  spaced  secondary  echoes, 
it  is  very  difficult  to  obtain  accurate  short  distance  measurements  for  acoustic  proxim- 
ity ranging  systems  using  traditional  matched  filter  based  methods.  In  this  chapter, 
a computationally  efficient  time  delay  estimation  algorithm,  referred  to  as  Multi- 
PEARS,  is  presented  for  the  joint  proximity  ranging  and  secondary  echo  mitigation. 
Numerical  results  demonstrate  that  Multi-PEARS  can  deal  with  secondary  echoes 
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Figure  6.5:  Recored  waveforms  in  the  real  experiment,  (a)  Received  signal  at  level 
3.  (b)  Received  signal  at  level  5.  (c)  Received  signal  at  level  15.  (d)  Received  signal 
at  level  19. 


effectively,  and  the  RMSEs  of  the  estimates  obtained  by  using  Multi-PEARS  can  ap- 
proach the  corresponding  CRBs  as  the  SNR  increases.  Experimental  results  obtained 
by  using  Panasonic  ceramic  transducers  show  that  Multi-PEARS  can  also  perform 
very  well  in  a practical  environment. 
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Sample  Points  Sample  Points 

(c)  (f) 

Figure  6.6:  Application  of  Multi-PEARS  to  a recorded  signal  at  level  19.  (a)  The 
original  signal,  (b)-(e)  The  reconstructed  signals  by  using  Multi-PEARS  estimates 
for  the  echoes  1 to  4,  respectively,  (d)  The  reconstructed  signal  by  superposition  of 
the  signals  in  (b)-(e). 
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Figure  6.7:  Stationary  level  measurement  results  obtained  by  applying  the  Multi- 
PE  ARS  algorithm  to  the  measured  signal  with  the  different  choices  of  L.  The  ranging 
distance  covers  from  29  mm  to  39  mm.  The  distance  between  two  consecutive  levels 
is  0.5  mm.  (a)  Measurement  results  vs.  the  true  distance  obtained  by  traverse,  (b) 
Measurement  errors. 


CHAPTER  7 

EXPERIMENTAL  RESULTS  IN  AIR- WATER  TUNNEL 

7.1  Introduction 

In  Chapters  5 and  6,  the  experiments  have  been  conducted  to  verify  the  per- 
formances of  the  proposed  distance  estimation  techniques.  In  those  experiments, 
the  smooth  aluminum  plate  is  used  to  simulate  the  sound-hard  boundary.  In  this 
chapter,  the  air-water  interfaces  in  the  tunnel,  which  more  resemble  the  real-world 
supercavitating  vehicle  measurement  environment,  are  measured  by  our  novel  prox- 
imity ranging  methods. 

The  photograph  of  the  experimental  setup  is  shown  in  Figure  7.1.  The  Pana- 
sonic ceramic  transducers  are  mounted  above  the  air- water  interface  and  their  geomet- 
ric centers  are  17  mm  apart.  The  same  excitation  and  transmitted  waveforms  adopted 
in  Chapter  6 (see  Section  6.4.1  for  details)  are  used  in  this  experiment.  The  excita- 
tion and  received  waveforms  are  recorded  by  the  VXI  system.  The  measurement  data 
are  transferred  to  a PC  and  processed  by  the  distance  estimation  algorithm.  Due  to 
the  fact  that  the  strong  secondary  echoes  can  be  induced  in  this  experimental  setup, 
the  Multi-PEARS  algorithm,  which  can  mitigate  the  secondary  echoes  effectively,  is 
employed  here  to  perform  the  distance  estimation.  In  the  experiment,  three  different 
air- water  interfaces  are  considered.  First,  the  stationary  interface  is  measured.  The 
results  are  used  to  show  the  measurement  accuracy  and  verify  the  assumption  of  the 
sound-hard  reflection.  Then,  the  nonstationary  air-water  interfaces  at  different  water 
flow  speeds  are  measured.  Finally,  the  experiment  is  performed  on  the  nonstationary 
interface  with  the  additional  periodic  disturbance  being  applied. 
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Figure  7.1:  Experimental  setup  of  the  acoustic  ranging  system  for  measuring  the 
air-water  interface  in  the  tunnel. 


7.2  Stationary  Air- Water  Interface 

In  this  case,  the  flow  speed  in  the  tunnel  is  zero  and  the  air-water  interface 
is  stationary  and  smooth.  The  sensors  are  fixed  on  the  traverse  which  is  vertically 
moving  with  a constant  speed  of  2 cm/s.  The  excitation  and  received  waveforms  are 
recorded  by  the  VXI  system  when  the  traverse  is  moving. 

The  distance  measurement  results  are  shown  in  Figure  7.2(a).  Two  methods 
are  tested  here.  One  is  the  Multi-PEARS  algorithm,  which  is  based  on  the  assump- 
tion that  the  amplitude  of  the  first  echo  is  nonnegative  and  is  denoted  by  dots,  and  the 
other  is  the  simplified  Multi-PEARS  algorithm,  which  is  based  on  the  complex-valued 
amplitude  assumption  and  is  denoted  by  circles.  (The  simplified  Multi-PEARS  algo- 
rithm uses  only  the  initialization  stage  of  the  Multi-PEARS  algorithm.  See  Section 
6.3.1.)  The  measurement  results  are  compared  to  the  true  distance  measured  by  the 
traverse,  and  the  errors  are  shown  in  Figure  7.2(b).  The  RMSEs  for  our  measure- 
ments are  0.045  mm  for  the  nonnegative  case  and  0.77  mm  for  the  complex-valued 
case  with  the  entire  measurement  range  covering  from  3 to  11  cm.  It  is  worth  not- 
ing that  although  the  measurement  error  of  the  simplified  Multi-PEARS  algorithm 
is  larger  than  that  of  the  Multi-PEARS  algorithm,  the  former  is  computationally 
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simpler  than  the  latter  since  there  is  no  need  for  the  iteration.  In  the  experiment, 
the  average  processing  time  for  one  estimate  using  Matlab  code  on  the  PC  (Pen- 
tium 4,  1.38  GHz)  is  about  0.63  s for  Multi-PEARS  and  0.36  s for  the  simplified 
Multi-PEARS.  Moreover,  in  the  scenarios  of  the  rough  surface  measurements  (i.e., 
when  the  irregularity  height  of  the  surface  is  on  the  order  of  the  wavelength),  the 
nonnegative  amplitude  assumption  is  no  longer  valid  and  only  the  complex-valued 
amplitude  assumption  based  method  can  be  used  [10,  112].  Therefore,  for  the  spec- 
ular sound-hard  surface  with  a high  ranging  accuracy  requirement,  Multi-PEARS 
should  be  used.  Otherwise,  the  simplified  Multi-PEARS  is  preferred. 

In  this  experiment,  the  attenuation  of  the  amplitude  of  the  first  echo  as  a 
function  of  the  distance  is  also  studied.  The  amplitude  measurement  of  the  first 
echo  is  estimated  by  the  simplified  Multi-PEARS  algorithm,  in  which  the  amplitude 
is  calculated  by  comparing  the  first  echo  with  the  reference  signal  whose  energy  is 
normalized.  The  absolute  value  of  the  amplitude  measurement  results  for  the  first 
echo  as  a function  of  the  distance  are  shown  in  Figure  7.2(c)  by  circles.  (The  results 
calculated  by  Multi-PEARS  are  very  similar  to  those  of  the  simplified  Multi-PEARS 
and  are  omitted  here.) 

The  relationship  between  the  echo  amplitude  and  the  distance  can  be  modeled 

as  [8] 

p(d)  = e~&^,  (7.1) 

d2 

where  the  term  e~ ^ is  used  to  model  the  directivity  of  the  acoustic  sensors  and  9 is 
defined  as  9 — tan_1(r/2d)  with  r being  the  distance  between  the  geometric  centers 
of  the  transmitter  and  the  receiver  and  d being  the  distance  between  sensors  surface 
to  the  air-water  interface  (see  Figure  1.1(b));  l/R  counts  for  the  geometry  dissipation 
of  the  acoustic  wave  propagating  in  the  free  space  with  R = 2i/r2/4  + d2;  a and  A0 
are  constants.  (A0  can  be  regarded  as  the  amplitude  measurement  when  R - 1 cm 
and  9 = 0°.)  Note  that  both  9 and  R are  functions  of  d.  Based  on  this  model, 


98 


12  - 

11  - 

_10  - 
E 

r 9- 

c 

© 

E 8- 

3 

</)  7 - 
<0  ' 

© 

l- 

2 5- 


,?»■ 


«?•* 


«*•*  Multi-PEARS  (Nonnegative) 

Simplified  Multi-PEARS  (Complex-Valued) 


_l | | | | L. 


1.5  2 2.5 

Time  (Second) 


(a) 


2.5  - 
2 - 

1.5  - 
1 - 


E 
iu 

c 0.5  - 
© 

E 

© Oh 


2 -0. 

2 


S 

</)  _1  c 

Q 


Multi-PEARS  (Nonnegative) 

Simplified  Multi-PEARS  (Complex-Valued) 


1.5  2 2.5 

Time  (Second) 


(b) 


(c) 


Figure  7.2:  Stationary  air-water  interface  measurements  in  the  tunnel.  The  traverse  is 
vertically  moving  with  a constant  speed  of  2 cm/s.  (a)  Distance  measurement  results 
based  on  different  assumptions  on  the  amplitude  of  the  first  echo,  (b)  Distance 
measurement  errors,  (c)  Amplitude  measurement  results. 


the  NLS  fitting  method  is  used  to  fit  the  amplitude  measurements,  and  we  obtain 
d = 16°  and  A0  = 6.3.  The  fitting  curve  is  shown  in  Figure  7.2(c)  by  the  solid  line. 

7.3  Nonstationary  Air- Water  Interface 

In  this  case,  the  height  of  the  sensors  with  respect  to  the  air-water  interface 
is  fixed  at  about  4 cm,  and  the  water  flow  speed  is  changed  by  increasing  the  control 
frequency  from  the  control  panel.  There  is  a linear  relationship  between  the  flow 
speed  and  the  control  frequency  as  shown  in  Figure  7.3,  where  the  circles  denote  the 
real  water  flow  speed  measurement  results  and  the  solid  line  is  obtained  by  the  least 
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Figure  7.3:  The  water  flow  speed  in  the  tunnel  as  the  function  of  the  frequency  set  by 
the  control  panel.  The  circles  denote  the  experimental  measurement  results  and  the 
solid  line  is  the  least  squares  fitting  result  based  on  the  experimental  measurements. 


squares  fitting  based  on  the  real  measurement  data.  The  linear  relationship  between 
the  water  flow  speed  and  the  control  frequency  can  be  expressed  as 


v = af, 


(7.2) 


where  a = 0.0143  m • s_1  • Hz-1  in  the  experiment. 

Then  we  measure  the  air-water  interface  by  acoustic  sensors  when  the  control 
frequencies  are  changed  from  0 Hz  to  60  Hz  with  an  increment  of  5 Hz.  The  distance 
and  the  amplitude  measurement  results  for  0 Hz  (0  m/s),  20  Hz  (0.286  m/s),  40  Hz 
(0.572  m/s),  and  60  Hz  (0.858  m/s)  are  shown  in  Figure  7.4.  The  entire  observation 
time  is  2 s,  in  which  382  estimates  for  both  the  distance  and  the  amplitude  can  be 
obtained.  All  of  the  measurement  results  shown  here  are  obtained  by  the  simplified 
Multi-PEARS  algorithm  due  to  the  wavy  air-water  interfaces.  Note  from  the  distance 
measurements  in  the  left  column  of  Figure  7.4  that  even  though  the  position  of 
the  sensors  is  fixed,  the  means  of  the  distance  measurements  are  slightly  different 
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for  the  different  water  flow  speeds.  This  is  because  that  as  the  water  flow  speed 
changes,  the  entire  air- water  interface  profile  in  the  tunnel  changes  as  well.  Thus, 
the  distance  between  the  sensors  and  the  local  air-water  interface  is  slightly  different 
at  the  different  water  flow  speed. 

According  to  (7.1),  the  amplitude  of  the  echo  is  a function  of  the  distance. 
Hence,  even  for  the  same  air-water  interface,  the  amplitude  measurements  can  be 
different  due  to  the  distance  difference  between  two  observations.  However,  such  a 
difference  in  the  amplitude  measurements  can  be  readily  compensated  based  on  (7.1) 
by  using  the  known  distance  measurements.  The  right  column  in  Figure  7.4  shows  the 
amplitude  measurements  after  the  distance  compensation  with  a reference  distance 
being  d0  — 4 cm. 

In  this  experiment,  it  is  noted  that  the  amplitude  measurements  are  generally 
more  sensitive  to  the  water  flow  speed  than  the  distance  measurements.  This  obser- 
vation inspires  us  to  propose  that  the  change  of  the  amplitude  measurements  could 
be  used  to  reflect  the  water  flow  speed  in  the  tunnel. 

Define  the  change  of  the  amplitude  measurements  as 


pensation,  and  N is  the  number  of  the  measurements.  The  change  of  the  amplitude 
measurements  as  a function  of  the  water  flow  speed  is  shown  in  Figure  7.5.  From 
this  figure,  one  can  clearly  see  the  trend  that  a increases  as  the  water  flow  speed 
increases. 

7.4  Nonstationary  Air- Water  Interface  with  Periodic  Disturbance 


(7.3) 


where  /3(n),  n = 1,2,  ...,N  are  the  amplitude  measurements  after  the  distance  com- 


Based on  the  same  experimental  setup  used  in  Section  7.3,  the  traverse  is 
further  used  to  generate  the  additional  water  wave  by  moving  a small  aluminum 
block  in  and  out  of  the  air- water  surface  continuously  with  a period  of  1.6  s.  Thus, 
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Time  (Second)  Time  (Second) 

Figure  7.4:  Air-water  interface  measurement  results  in  the  tunnel  with  different  water 
flow  speeds.  Left  column:  Distance  measurement  results.  Right  column:  Amplitude 
measurement  results. 

in  this  setup,  there  are  two  factors  contributing  to  the  irregularity  of  the  air-water 
interface  (i.e.,  the  constant  motion  of  the  water  flow  and  the  periodic  disturbance 
generated  by  the  traverse). 

The  distance  and  the  amplitude  measurements  when  the  water  flow  speed 
v - 0.43  m/s  (i.e.,  the  control  frequency  / = 30  Hz)  are  shown  in  Figures  7.6(a)  and 
7.6(b),  respectively.  The  observation  time  is  16  s.  Note  that  the  periodic  behavior  due 
to  the  disturbance  can  be  clearly  observed  from  both  the  distance  and  the  amplitude 
measurements.  The  discrete  Fourier  transform  of  the  distance  and  the  amplitude 
measurements  (with  the  DC  component  removed)  are  shown  in  Figures  7.6(c)  and 
7.6(d),  respectively.  The  dash-dot  lines  denote  the  true  location  of  the  disturbance 
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Figure  7.5:  The  change  of  the  amplitude  measurements  as  the  function  of  the  water 
flow  speed. 


frequency.  The  measurement  results  for  the  water  flow  speed  v — 0.57  m/s  (i.e.,  the 
control  frequency  / = 40  Hz)  are  shown  in  Figure  7.7. 

7.5  Summary 

In  this  chapter,  the  different  air- water  interfaces  in  the  tunnel  at  different  water 
flow  conditions  are  measured  by  the  acoustic  proximity  ranging  system.  First,  the 
stationary  air-water  interface  is  measured.  It  is  shown  that  the  distance  measurement 
accuracy  can  achieve  0.045  mm  for  the  Multi-PEARS  algorithm  (i.e.,  the  nonnegative 
amplitude  assumption  based  method)  and  0.7  mm  for  the  simplified  Multi-PEARS 
algorithm  (i.e.,  the  complex- valued  amplitude  assumption  based  method).  Then  the 
nonstationary  air- water  interfaces  with  the  different  water  flow  speeds  are  measured. 
The  change  of  the  amplitude  measurements  is  also  introduced  to  indicate  the  flow 
speed.  Finally,  the  nonstationary  interfaces  with  the  additional  periodic  disturbance 
are  measured.  The  experimental  results  show  that  the  periodic  disturbance  can  be 
correctly  measured  by  the  proximity  ranging  system. 
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Figure  7.6:  Air-water  interface  measurement  results  when  the  water  flow  speed 
v = 0.43  m/s  (i.e.,  / = 30  Hz)  and  the  disturbance  period  is  1.6  s.  (a)  Distance 
measurement  results,  (b)  Amplitude  measurement  results,  (c)  Discrete  Fourier  trans- 
form of  the  distance  measurements  with  the  DC  component  removed.  (The  dash-dot 
line  denotes  the  true  frequency  location  of  the  disturbance.)  (d)  Discrete  Fourier 
transform  of  the  amplitude  measurements  with  the  DC  component  removed. 
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Figure  7.7:  Air- water  interface  measurement  when  the  water  flow  speed  v = 0.57 
m/s  (i.e.,  / = 40  Hz)  and  the  disturbance  period  is  1.6  s.  (a)  Distance  measurement 
results,  (b)  Amplitude  measurement  results,  (c)  Discrete  Fourier  transform  of  the 
distance  measurements  with  the  DC  component  removed.  (The  dash-dot  line  denotes 
the  true  frequency  location  of  the  disturbance.)  (d)  Discrete  Fourier  transform  of  the 
amplitude  measurements  with  the  DC  component  removed. 


CHAPTER  8 
CONCLUSIONS 

8.1  Summary  Remarks 

In  this  work,  we  have  studied  using  acoustic  proximity  ranging  techniques 
for  monitoring  the  cavity  thickness  of  a supercavitating  vehicle.  Specifically,  the 
phase-shift  approach,  the  multi- frequency  technique,  the  PEARS  scheme,  and  the 
Multi-PEARS  algorithm  are  proposed  for  such  a purpose.  Both  theoretical  analyses 
and  experiments  conducted  by  the  commercially  available  acoustic  transducers  are 
used  to  verify  the  ranging  performance  of  the  new  methods. 

In  Chapter  3,  we  begin  this  work  by  the  theoretical  study  of  the  phase-shift 
approach,  which  is  one  of  the  simplest  ranging  methods.  For  the  phase-shift  approach, 
we  assume  that  only  two  frequencies  are  used  and  the  measurements  are  corrupted 
by  colored  Gaussian  noise.  By  taking  the  a priori  knowledge  of  the  acoustically 
hard  reflection  into  account,  a new  time  delay  estimation  algorithm  based  on  the 
ML  theory  is  derived.  Numerical  examples  show  that  our  new  method  outperforms 
the  traditional  method  in  terms  of  both  the  estimation  accuracy  and  the  robustness 
against  data  model  mismatch. 

A straightforward  extension  of  the  phase-shift  approach  is  addressed  in  Chap- 
ter 4,  where  an  arbitrary  number  of  frequencies  is  considered.  In  this  chapter,  first, 
a new  multi-frequency  CW  waveform  is  designed  to  simplify  the  procedures  of  the 
phase-shift  measurement.  Then,  a novel  time  delay  estimator  based  on  the  NLS  fit- 
ting criterion  is  derived.  To  minimize  the  highly  oscillatory  cost  function,  an  efficient 
two-stage  estimation  algorithm  is  proposed.  It  is  shown  by  numerical  examples  that 
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for  a fixed  frequency  interval  and  a fixed  SNR,  the  more  frequencies  we  use,  the  lower 
the  SNR  threshold  and  the  higher  the  estimation  accuracy  can  be  obtained. 

Inspired  by  the  multi-frequency  technique  studied  in  Chapter  4,  the  so-called 
PEARS  scheme  is  proposed  in  Chapter  5.  PEARS  is  novel  in  that  it  is  applicable 
to  arbitrary  transmitted  waveforms  as  long  as  the  waveform  is  periodic.  Therefore, 
PEARS  can  provide  more  flexibility  in  designing  a practical  proximity  ranging  system. 
Experimental  results  obtained  by  the  commercially  available  ultrasonic  transducers 
demonstrate  the  good  performance  of  PEARS.  However,  in  the  experiments,  we  also 
notice  that  the  performance  of  our  new  methods  may  degrade  due  to  the  presence  of 
the  secondary  echoes.  The  problem  noted  in  the  experiments  leads  us  to  the  topic  of 
Chapter  6. 

To  deal  with  the  secondary  echoes,  in  Chapter  6,  a new  data  model  and 
the  corresponding  time  delay  estimation  algorithm,  referred  to  as  Multi-PEARS,  is 
presented  for  the  joint  proximity  ranging  and  secondary  echo  mitigation.  The  key  idea 
of  the  Multi-PEARS  algorithm  is  to  obtain  a good  initial  time  delay  estimate  using 
the  FFT  based  fast  search  technique,  and  then,  refine  the  initial  estimate  by  using 
the  former  developed  PEARS  algorithm.  Both  numerical  and  experimental  results 
demonstrate  that  Multi-PEARS  can  provide  very  accurate  distance  measurements 
even  in  the  presence  of  the  strong  secondary  echoes. 

Since  this  work  is  highly  application-oriented,  the  ranging  experiments  are 
further  conducted  in  the  air-water  tunnel,  and  the  results  are  illustrated  in  Chapter  7. 
Three  different  air- water  interfaces  (i.e. , stationary,  nonstationary,  and  nonstationary 
with  the  periodic  disturbance)  are  measured.  Experimental  results  show  that  our 
acoustic  proximity  ranging  technique  can  provide  the  accurate  distance  information 
between  sensors  and  air-water  interface  at  different  flow  conditions  in  the  tunnel. 

A brief  summary  of  the  acoustic  proximity  ranging  techniques  proposed  in 
this  dissertation  are  illustrated  in  Table  8.1. 
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Table  8.1:  Summary  of  the  different  acoustic  proximity  ranging  techniques  proposed 
in  this  dissertation. 


Phase-shift 

Multi-frequency 

PEARS 

Multi-PEARS 

Waveform 

CW;  Two 
frequencies. 

CW;  Arbitrary 
number  of 
frequencies; 
Periodic. 

Arbitrary; 

Periodic. 

Arbitrary; 

Periodic. 

Usage  of  Transducer 
Bandwidth 

Low 

High 

High 

High 

Noise  Assumption 

Colored 

White 

White 

White 

Pre-Processing 

IQ  or  FFT 

FFT 

FFT 

FFT 

Time  Delay 
Estimate 

Initial 

Closed-form 

Closed-form 

FFT-based 

search 

FFT-based 

search 

Refined 

Closed- form 

Closed-form 

Closed-form 

Closed-form  + 
Iteration 

Robustness  Against 
Secondary  Echoes 

No 

No 

No 

Yes 

8.2  Future  Work 

The  future  work  can  be  divided  into  two  aspects  of  the  theoretical  study  for  the 
proximity  ranging  algorithm  and  the  implementation  of  the  proposed  techniques.  In 
the  theoretical  study  aspect,  two  interesting  areas  closely  related  to  the  current  work 
can  be  investigated.  The  first  case  would  be  the  analysis  of  the  multi-frequency  tech- 
nique or  the  PEARS  scheme  in  the  colored  noise.  Unlike  the  phase-shift  approach, 
there  may  be  no  explicit  expression  to  the  time  delay  estimate.  The  search  based  op- 
timization methods  will  be  considered  to  maximize  the  corresponding  ML  functions. 
Another  interesting  study  would  be  the  extension  of  the  Multi-PEARS  algorithm  by 
assuming  a constraint  on  the  amplitude  for  each  echo.  That  is,  the  amplitude  of  the 
echo  should  decease  as  the  time  delay  of  the  echo  increases.  This  constraint  is  fairly 
reasonable  due  to  the  geometrical  dissipation  of  the  sound  wave  propagating  in  the 
media.  By  using  this  constraint,  the  time  delay  estimation  accuracy  is  expected  to 
improve  for  a low  to  moderate  SNR. 
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In  the  implementation  aspect,  first,  the  MEMS  acoustic  proximity  transducers 
will  be  tested.  It  is  important  to  find  out  the  key  parameters  of  the  MEMS  trans- 
ducer such  as  SNR  range,  transmitting  power,  sensitivity,  and  directivity,  etc.  With 
these  parameters  in  hand,  the  best  suited  proximity  ranging  approach  can  be  chosen 
from  our  former  developed  methods.  Then  the  corresponding  distance  estimation 
algorithm  will  be  programmed  into  the  high  speed  DSP  chips.  Finally,  the  prototype 
of  the  MEMS  based  real-time  acoustic  ranging  system  will  be  built  and  tested. 


APPENDIX  A 

THE  DERIVATION  OF  CRB  FOR  PHASE-SHIFT  APPROACH 

We  outline  below  the  derivation  of  CRBs  for  the  parameter  estimates  of  the 
following  data  model: 


x(n)  = (3b  + e(n),  n = 1, 2, N,  (A.l) 

where  b = [ ej27r^lT  ej27rAT 

In  (A.l)  the  additive  noise  vectors  e(n),  n = 1,2, ...,  N are  assumed  to  be  inde- 
pendent zero-mean  Gaussian  random  vectors  with  an  unknown  covariance  matrix  Q. 
The  data  samples  x(n),  n — 1,2, ...,  N are  independent  and  the  unknown  parameters 
in  the  estimation  problem  are  the  real  and  imaginary  parts  of  (3  (for  complex-valued 
/?)  or  simply  the  amplitude  of  /3  (for  real-valued  /3),  the  time  delay  r and  the  unknown 
elements  of  Q.  Let  77  be  a vector  containing  those  parameters.  The  CRB  inequality 
states  that  provided  the  estimates  are  unbiased,  the  covariance  matrix  of  the  esti- 
mated parameter  vectors  77  is  lower  bounded  by  the  inverse  of  the  Fisher  information 
matrix  FIM-1.  The  FIM  can  be  computed  by  the  extended  Slepian-Bang’s  formula 
[99]: 

= WtrfQ-'QjQ-'Q')  + 2AfRe[(/Jb)'"Q-1  (/3b);],  (A.2) 

(here  X-  denotes  the  derivative  of  X with  respect  to  the  Ah  unknown  parameter). 
Note  that  the  FIM  is  a block  diagonal  matrix  since  Q does  not  depend  on  the  param- 
eters in  /3b (r),  and  f3b(r)  does  not  depend  on  the  elements  in  Q.  Hence,  the  CRB 
of  the  estimates  of  the  delay  and  the  amplitude  can  be  determined  from  the  second 
term  of  the  right  side  of  (A.2). 
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no 


We  first  derive  the  CRBs  for  the  case  of  a complex-valued  /3.  Let 


rjsc  = [ Re(/3)  Im(/?)  r ], 


(A.3) 


where  Im(x)  denotes  the  imaginary  part  of  x,  be  a vector  of  the  signal  parameters 
(i.e.,  the  parameters  in  the  upper  left  block  of  the  FIM).  Define 


and  let 


Then 


II 

CD 

II 

cr 

(A.4) 

dlm{p)  3 ’ 

(A.5) 

' dV  -o  o r/ie^Ar- 

dr  32n^  [/2eJ'2,r/^J  ’ 

(A.6) 

F = [ A*i  A 4 Ms  ]• 

(A-7) 

CRBfoJ  = [2AfRe(FwQ-1F)]-1, 

(A.8) 

which  is  easily  evaluated  numerically.  If  a2  — of  = o\  and  p — 0,  simplifying  (A. 8) 
yields 

CRB„  = CRB^,  + CRBIm((J)  = ^(1  + 2 /j+jL),  (A.9) 

and 


CRBt  = 


(A-10) 


4AT|/J|  V(h-h)2 

Clearly,  for  maximal  accuracy,  /2  — f\  should  be  as  large  as  possible.  Note  however 
that  f2  — f\  < (cf.  the  remark  before  (3.42)). 

Second,  we  derive  the  CRBs  for  real-valued  /3.  Let 


Define 


%,  = [ P r]. 


(A.ll) 


(A.12) 


Ill 


(A.13) 


Let 


F = [ iA  Ih  ]• 


(A.14) 


(A-15) 


(A. 16) 


and 


CRB'  sWi'l/ft/l)' 


(A.  17) 


In  general,  the  frequency  difference  is  determined  by  the  maximum  distance 
we  want  to  measure.  Therefore,  once  |/x  — /2|  is  fixed,  we  can  choose  large  f\  and  /2 
to  decrease  the  real  CRB  of  r.  Note  however  that  the  lower  real  CRB  we  design,  the 
higher  SNR  is  required  to  approach  it. 

Note  that  in  the  case  when  Q = a2 1,  according  to  (A. 10)  and  (A. 17),  the  ratio 
of  the  CRBs  of  r for  complex-valued  /3  and  for  real-valued  /?  is 


This  confirms  the  heuristic  reasoning  in  Section  3.4  which  concludes  that  the  cost 
function  with  the  sharpest  peak  should  give  the  most  accurate  estimate.  (Note  that  in 
the  numerical  examples,  for  convenience,  the  square  roots  of  the  CRBs  are  illustrated 
in  the  corresponding  figures.) 


(A.  18) 


This  ratio  is  directly  related  to  the  ratio  of  the  oscillation  frequencies  of  C\(r)  and 
c2(t),  respectively,  and  therefore  to  the  sharpness  of  the  peaks  of  C\{r)  and  c2(r). 


APPENDIX  B 

THE  DERIVATION  OF  CRB  FOR  MULTI-FREQUENCY  TECHNIQUE 

We  outline  below  the  derivation  of  CRBs  for  the  parameter  estimates  of  the 
data  model  in  (4.10).  Define 

y = [y(0)  y{  1)  ...  y(N  - 1)  ]T,  (B.l) 

e = [ e(0)  e(l)  ...  e(N  - 1)  ]T,  (B.2) 

and 

g(ir)  = [p(0)  P(l)  - s(Y-l)]T,  (B.3) 

where 

K- 1 

9(n)  = ^ \HUk)\cos{2nfkn-2TTfkT+^k+avg[H(fk)]},  n = 0, 1, N-l.  (B.4) 

k=0 

Then  (4.10)  can  be  expressed  as 

y = /3g(r)  + e,  (B.5) 

where  the  additive  noise  vector  e is  assumed  to  be  a zero-mean  real-valued  Gaussian 
random  vector  with  a covariance  matrix  Q = a2 1 with  I being  the  identity  matrix. 
The  unknown  parameters  in  the  estimation  problem  are  the  real-valued  /3,  the  time 
delay  r,  and  the  unknown  variance  a2.  Let  r]  be  a vector  containing  these  parameters. 
The  CRB  of  r]  is  the  inverse  of  the  Fisher  information  matrix  FIM,  where  the  (zj)th 
element  of  FIM  can  be  computed  by  the  extended  Slepian-Bang’s  formula  [99]: 

{FIM},;  = itr(Q-'Q;Q-‘Q')  + [(£tg);rQ-‘(/3g)' ],  (B.6) 

and  here  Xj  denotes  the  derivative  of  X with  respect  to  the  zth  unknown  parameter. 
Note  that  FIM  is  a block  diagonal  matrix  since  Q does  not  depend  on  the  parameters 
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in  /3g (r),  and  /3g (r)  does  not  depend  on  the  elements  in  Q.  Hence,  the  CRB  of  the 
time  delay  and  the  amplitude  can  be  determined  from  the  second  term  of  the  right 
hand  side  of  (B.6). 

Let 

V , = [ r & ]T,  (B.7) 

be  the  vector  containing  the  signal  parameters.  Define 


AH 


AH 


djg 

dp 

dpg 

dr 


= g, 


= , 


where 


s'  = = r M21  Mil  gg(Jv— i)  it 

6 1 0T  9r  dr  j ’ 


with 


dr 


K- 1 


= 27Tfk\H(fk)\ sm{2nfkn  - 2tt fkr  + ipk  + arg [#(/*)]}. 


k=0 


Let 


p = [ /*i 


(B.8) 

(B.9) 

(B.10) 

(B.ll) 

(B-12) 


Then 

CRBfa.)  = (FhQ-1F)-1,  (B.13) 


which  is  easily  evaluated  numerically. 


APPENDIX  C 

THE  DERIVATION  OF  CRB  FOR  PEARS  SCHEME 


We  sketch  below  the  derivation  of  the  CRBs  for  the  real-valued  parameters 
of  the  data  model  in  (5.3).  Due  to  the  conjugate  symmetric  property  of  DFTs, 


{Y{k)}k--N/2  can  be  expressed  in  terms  of  {Y(k)}kl_N/2 
complex-valued  and  {Y(-N/2),  Y(0)}  being  real-valued. 

with  {Y{k)}ki_ 
Let 

-n/  2+1  being 

Yc  = [ Y(—N/2  + 1)  Y{-N/2  + 2)  ••• 

Y{- 1)  F, 

(C.l) 

Yr  = [ Y(-N/2)  F(0)  ]T, 

(C.2) 

Ec  = [ E(—N/2  + 1)  E{-N/ 2 + 2)  ••• 

E(- 1)  ]T, 

(0.3) 

Er  = [ E(—N/2)  E{ 0)f, 

(0.4) 

Sc  = [ S(—N/2  + 1)  S{-N/2  + 2)  • • • 

S{- 1)  F, 

(C.5) 

Sr  = [ S(—N/2)  5(0)  ]r, 

(C.6) 

bc(r)  = [ ••• 

e-^(-p  f, 

(C.7) 

and 

br(r)  = [ 1 ]T- 

(C.8) 

We  have 

Yc  = + Ec, 

(C.9) 

and 

Yr  = fir  + Er, 

(C.10) 

where 

H-c  — /^Scbc, 

(C.ll) 

114 


115 


and 

llr  = PSrbr.  (C.12) 

Assume  that  the  additive  noise  {e(nT5)}^T01  is  a real- valued  zero-mean  white 
Gaussian  random  process  with  variance  a2.  Denote 

0 = [ P t a2  ]T.  (C.13) 

Since  the  DFT  matrix  is  a unitary  operator,  the  joint  probability  density  function 
for  {Yc,  Yr}  has  the  form: 

p({Yc,Yr}|0)  = - ^(Yc  - mJH(Yc  - Mc) 

~2j2(Yr  _ /Or(Yr  -/*!■)}•  (C.14) 

Based  on  the  extended  Slepian-Bang’s  formula,  the  CRB  matrix  can  be  com- 
puted as 

[CRB-'Wli,  = ^.r[Q,lQ,;qr-1Q,']  + [pr;7Q,-1Pr'] 

trlQe-'Qc'Qr'Qc;]  + 2Re[/xc;',Qc-Vci],  (C15> 

where  Qr  = o2lr  with  Ir  being  2x2  identity  matrix,  Qc  = <t2Ic  with  Ic  being 
(N/2  - 1)  x (N/2  - 1)  identity  matrix,  and  X-  denotes  the  derivative  of  X with 
respect  to  the  zth  unknown  parameter.  Since  Qr  and  Qc  do  not  depend  on  the 
parameters  in  fir  and  nc,  and  /zr  and  /ic  do  not  depend  on  the  elements  in  Qr  and 
Qc,  it  can  be  shown  that  the  matrix  CRB(0)  is  block  diagonal  with  its  last  row  and 
last  column  being  zero  except  for  the  last  diagonal  element.  Let  the  signal  parameter 
vector  77  be  denoted  as 

V = [ P TV-  (C.16) 

Then 

[CRB-1  (»7)]*j  = [/xr;TQr-Vri]  + 2Re[/ic;HQc"Vci-] 

= 2Re[(i'HCrVj]. 


(C.17) 
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where 


li  = ^Wb(r), 


and  Q = a2 1 with  I being  the  identity  matrix  of  dimension  N/2  + 1. 


Define 


where 


dfx 

dp 

dfx 

dr 


--  Wb(r), 
^Wb'(r), 


27T 


b'M  = -iw l (-f  + 


Of 


and 


(C.18) 

(C.19) 

(C.20) 

(C.21) 

(C.22) 


Then 


CRB  (r/)  = ytRe^F)]-1. 


(C.23) 


APPENDIX  D 

THE  DERIVATION  OF  CRB  FOR  MULTI-PEARS  ALGORITHM 


We  sketch  below  the  derivation  of  the  CRBs  for  the  data  model  in  (6.3) 
where  {a/}^  are  assumed  to  be  real-valued.  Due  to  the  conjugate  symmetry  prop- 
erty of  DFTs,  {Y{m)}^ZlN/ 2 can  be  expressed  in  terms  of  {Y(m)}Qm=_N/2  with 
{^(m)}m1=-w/2+ 1 being  complex-valued  and  {Y(-N/2),Y(0)}  being  real-valued.  Let 


Yc  = [ Y(—N/2  + 1)  Y(—N/2  + 2)  •••  Y(-l)  ]T,  (D.l) 

Yr  = [ Y(—N/2)  Y(0)  ]T,  (D.2) 

Ec  = [ E(—N/2  + 1)  E(—N/2  + 2)  •••  £(-l)  ]T,  (D.3) 

Er  = [ E(-N/2)  E( 0)  ]T,  (D.4) 

Sc  = [ S(-JV/2  + l)  S(—N/2  + 2)  •••  S(- 1)  ]T,  (D.5) 

Sr  = [ S(—N/2)  5(0)  ]T,  (D.6) 

g(/r  + A,)  = [ e-^(iT+A')(-£)  e-i^(/r+A,)(-f+i)  ...  l ]T  (D.7) 

gc(/r  + A,)  = [ e-^(iT+Ai)(-£+1)  e-:,^(ZT+A,K-f+2)  ...  e~i wfr(iT+ am-i)  jt 

(D.8) 

gr(Ir  + Aj)  = [ e-J7%(<T+Ad(-f)  i ]T,  (D.9) 

G = [g(T)  g(2r  + A2)  •••  g(Lr  + AL)],  (D.10) 

Gc  = [gc(r)  gc(2r  4-  A2)  •••  gc(Lr  + AL)  ],  (D.ll) 

and 

Gr  = [gr(r)  gr(2r  + A2)  •••  gr(Z/r  + AL)  ].  (D.12) 
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We  have 

Yc  = nc  + EC|  (D.13) 

and 

Yr  — A*r  + Er,  (D.14) 

where 

AG  = ScGca,  (D.15) 

and 

AG  = SrGrc*,  (D.16) 

with 

oc  = [ ai  a2  aL  ]T.  (D.17) 

Assume  that  the  additive  noise  {e(nTs)}n=o  is  a real- valued  zero-mean  white 
Gaussian  random  process  with  variance  a2.  Denote 

0 = [ a i •••  aL  t A2  •••  Al  a2]T.  (D.18) 

Since  the  DFT  matrix  is  a unitary  operator,  the  joint  probability  density  function 
for  {Yc,  Yr}  has  the  form: 

P({Yc,Yr}|0)  = ^72aJVexp{  - ^(Y c - fic)H(Yc  - /xc) 

~2^(Yr~^)T(Yr-/ir)}-  (D-!9) 

Based  on  the  extended  Slepian-Bang’s  formula,  the  CRB  matrix  can  be  com- 
puted as 

[CRB-'Wfe  = 5tr[Qr-‘Qr'Qr-'Qr;]  + [A.r;TQrV,'l 

trlQe-'QwQ^Q,']  + 2Re[,ic;''QrlMc;],  (D.20) 

where  Qr  = cr2Ir  with  Ir  being  the  2 x 2 identity  matrix,  Qc  = a2 Ic  with  Ic  being  the 
(tV/2  — 1)  x (N/ 2 — 1)  identity  matrix,  and  Xj  denotes  the  derivative  of  X with  respect 
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to  the  ith  unknown  parameter.  Since  Qr  and  Qc  do  not  depend  on  the  parameters 
in  nr  and  fic , and  fir  and  nc  do  not  depend  on  the  elements  in  Qr  and  Qc,  it  can  be 
shown  that  the  matrix  CRB(0)  is  block  diagonal  with  its  last  row  and  last  column 
being  zero  except  for  the  last  diagonal  element.  Let  the  signal  parameter  vector  77 
be  denoted  as 

77  = [ ai  •••  aL  t A2  •••  A l]T-  (D.21) 


Let 


H = WGa,  (D.22) 

with  W being  defined  in  (7),  and  Q = a2 1 with  I being  the  identity  matrix  of 
dimension  N/2  + 1.  Then 


[CRB-1  (77)]^  = [Mr;TQr-1^]  + 2Re[/xc'//Qc-1^] 


vri  '‘vr  r~rj] 

fH^_  1 /• 

4jj 


= 2Re[/x'iWQ  V']. 


Define 


d/Jt 

dai 


= Wg(/r  + A,),  / = 1, 2, ...,  L, 

dfi  W3G 

a7  = wa7“' 


where 


with 


_ r dg(r)  dg(2r+A2) 
Qj-  1 dr  dr 


dg(Lr+AL)  1 
dr  b 


(D.23) 

(D.24) 

(D.25) 

(D.26) 


dg (It  + A;)  . 27 r/ 


dr 

and 

where 

dg  (It  + A;) 
dA  1 


(_|)e-i*(^A.K-¥)  (_|  + ...  o]T, 

ly  1 s 


dp  ^ Tdg(lr  + A/)  , 00  T 

= W — cti , l = 2,6, ...,  L, 


dA 


dAt 


(D.27) 


(D.28) 


-j^L[  (-^)e-i^(,T+^)(-f)  (_A  + i)e-*(^,)(-f +D  ...  0]r 

(D.29) 
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and 


F r djl  dll  dll  dll  dll  i 

l da\  dai  dr  dA^  dA £ 


(D.30) 


Then 


CRBfa)  = y[Re(F"F)]-1. 


(D.31) 
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